Commit 76e9e8ed authored by Michal Bursa's avatar Michal Bursa
Browse files

disk model for NT

parent 67ab5509
...@@ -3,15 +3,16 @@ ...@@ -3,15 +3,16 @@
CC = gcc CC = gcc
#SDLIB = /home/bursa/projects/slimdisk SIM5LIB = sim5lib
SIM5LIB = /home/bursa/projects/sim5/lib2
#CFLAGS = -I$(SDLIB) -I$(SIM5LIB) -L$(SDLIB) -Wall -O3 #SDLIB = /home/bursa/projects/slimdisk
#SIM5LIB = /home/bursa/projects/sim5/lib
#FLAGS = -I$(SDLIB) -I$(SIM5LIB) -L$(SDLIB) -L$(SIM5LIB) -Wall -O3
CFLAGS = -I$(SIM5LIB) -L$(SIM5LIB) -Wall -O3 -w CFLAGS = -I$(SIM5LIB) -L$(SIM5LIB) -Wall -O3 -w
LFLAGS = $(CFLAGS) -lm -ldl -lgsl -lgslcblas LFLAGS = $(CFLAGS) -lm -ldl -lgsl -lgslcblas
default: sd-radial-model default: disk-model
# compile object files as a "position independent code" # compile object files as a "position independent code"
...@@ -34,6 +35,8 @@ sd-radial-model: $(sd-radial-model-obj) ...@@ -34,6 +35,8 @@ sd-radial-model: $(sd-radial-model-obj)
sd-opacity-src = \ sd-opacity-src = \
$(SDLIB)/fluxmodel.c \ $(SDLIB)/fluxmodel.c \
$(SIM5LIB)/sim5kerr.c \ $(SIM5LIB)/sim5kerr.c \
......
...@@ -52,7 +52,7 @@ int main(int argc, char *argv[]) ...@@ -52,7 +52,7 @@ int main(int argc, char *argv[])
nRS = nrad*100; nRS = nrad*100;
RS = (radial_structure*) calloc(nRS, sizeof(radial_structure)); RS = (radial_structure*) calloc(nRS, sizeof(radial_structure));
disk_nt_setup(M, a, lum, alpha, 0); disk_nt_setup(M, a, lum, alpha, DISK_NT_OPTION_LUMINOSITY);
rmin = disk_nt_r_min(); rmin = disk_nt_r_min();
fprintf(stdout, "# radial disk structure model\n"); fprintf(stdout, "# radial disk structure model\n");
...@@ -66,6 +66,8 @@ int main(int argc, char *argv[]) ...@@ -66,6 +66,8 @@ int main(int argc, char *argv[])
fprintf(stdout, "# N_rad: %d\n", nrad); fprintf(stdout, "# N_rad: %d\n", nrad);
fprintf(stdout, "# model mdot: %e\n", disk_nt_mdot()); fprintf(stdout, "# model mdot: %e\n", disk_nt_mdot());
fprintf(stdout, "# model lumi: %e\n", disk_nt_lum()); fprintf(stdout, "# model lumi: %e\n", disk_nt_lum());
fprintf(stdout, "# format: ring-ID, radius [rg], temp [K], col.density [g/cm2], vert.grav. [?], vert.height [rg], spec.ang.mom. [], rad.velocity [c]\n");
// evaluate local model parameters at r=r0 // evaluate local model parameters at r=r0
int ir; int ir;
...@@ -109,35 +111,30 @@ int main(int argc, char *argv[]) ...@@ -109,35 +111,30 @@ int main(int argc, char *argv[])
RS[ir].Q = M*solar_mass*grav_const/pow(RS[ir].R*M*grav_radius,3.) * (sqr(u_f) + a*a*(u_t-1.))/RS[ir].R; RS[ir].Q = M*solar_mass*grav_const/pow(RS[ir].R*M*grav_radius,3.) * (sqr(u_f) + a*a*(u_t-1.))/RS[ir].R;
} }
fprintf(stderr, "L0: %e / %.3f\n", RS[nRS-1].totL, log(RS[nRS-1].totL)); double dL = (RS[nRS-1].totL-RS[1].totL);
fprintf(stderr, "L1: %e / %.3f\n", RS[1].totL, log(RS[1].totL));
double dL = (RS[nRS-2].totL-RS[1].totL);
fprintf(stderr, "Lmax=%.3f dL: %.3f\n", (RS[nRS-2].totL-RS[1].totL), dL);
int ID; int ID;
double prevL, prevR; double prevL, prevR;
do { do {
dL /= 1.1; dL /= 1.1;
ID = 0; ID = 0;
prevL = prevR = 0.0; prevL = prevR = 0.0;
for (ir=0; ir<nRS; ir++) { for (ir=0; ir<nRS; ir++) {
if ((ir==0) || (ir==nRS-1) || ((RS[ir].totL-prevL >= dL)&&(RS[ir].R-prevR>0.2)) || ((RS[ir].R-prevR)/RS[ir].R>0.5)) { if ((ir==0) || (ir==nRS-1) || ((RS[ir].totL-prevL >= dL)&&((RS[ir].R-prevR)/RS[ir].R>0.01)) || ((RS[ir].R-prevR)/RS[ir].R>0.5)) {
ID++; ID++;
prevL = RS[ir].totL; prevL = RS[ir].totL;
prevR = RS[ir].R; prevR = RS[ir].R;
} }
} }
fprintf(stderr, "dL= %e %d/%d\n", ID, nrad); } while (ID<=nrad);
} while (ID!=nrad);
ID = 0; ID = 0;
prevL = prevR = 0.0; prevL = prevR = 0.0;
for (ir=0; ir<nRS; ir++) { for (ir=0; ir<nRS; ir++) {
if ((ir==0) || (ir==nRS-1) || ((RS[ir].totL-prevL >= dL)&&(RS[ir].R-prevR>0.2)) || ((RS[ir].R-prevR)/RS[ir].R>0.3)) { if ((ir==0) || (ir==nRS-1) || ((RS[ir].totL-prevL >= dL)&&((RS[ir].R-prevR)/RS[ir].R>0.01)) || ((RS[ir].R-prevR)/RS[ir].R>0.5)) {
fprintf( fprintf(
stdout, "%04d/%04d %.2e %.4e %e %e %e %e %e %e\n", stdout, "%04d %.2e %.4e %e %e %e %e %e %e\n",
ID+1, ir, RS[ir].totL, RS[ir].R, RS[ir].T, RS[ir].S, RS[ir].Q, RS[ir].H, RS[ir].l, RS[ir].V ID+1, RS[ir].R, RS[ir].T, RS[ir].S, RS[ir].Q, RS[ir].H, RS[ir].l, RS[ir].V
); );
ID++; ID++;
prevL = RS[ir].totL; prevL = RS[ir].totL;
......
...@@ -39,7 +39,7 @@ int main(int argc, char *argv[]) ...@@ -39,7 +39,7 @@ int main(int argc, char *argv[])
int nrad = atoi(argv[5]); // number of rings int nrad = atoi(argv[5]); // number of rings
// initialize flux model // initialize flux model
flux_model_init("./flux-sd.so", M, a, lum, alpha, NULL, SD_INIT_OPT_LUMINOSITY); flux_model_init("./flux-nt.so", M, a, lum, alpha, NULL, SD_INIT_OPT_LUMINOSITY);
int ir; int ir;
for (ir=0; ir<nrad; ir++) { for (ir=0; ir<nrad; ir++) {
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment