Commit b29684e5 authored by Michal Bursa's avatar Michal Bursa
Browse files

added disk model

parent cbb77131
/data
/disk/*
!/disk/Makefile
!/disk/*.c
!/disk/*.h
/tlusty/*
!/tlusty/t205*.f
!/tlusty/data
......@@ -8,3 +13,5 @@
/plots/*
!/plots/*.asy
!/plots/*.py
/tmp
\ No newline at end of file
# $Id: Makefile.def,v 1.3 2005/07/18 12:58:34 bursa Exp $
CC = gcc
SDLIB = /home/bursa/projects/slimdisk
SIM5LIB = /home/bursa/sim5/lib
CFLAGS = -I$(SDLIB) -I$(SIM5LIB) -L$(SDLIB) -Wall -O3
LFLAGS = $(CFLAGS) -lm -ldl -lgsl -lgslcblas
default: sd-radial-model
# compile object files as a "position independent code"
%.o: %.c
$(CC) $(CFLAGS) -fpic -c -o $@ $<
clean:
rm -f *.o
sd-radial-model-src = \
$(SDLIB)/fluxmodel.c \
$(SIM5LIB)/sim5kerr.c \
sd-radial-model.c
sd-radial-model-obj = $(sd-radial-model-src:.c=.o)
sd-radial-model: $(sd-radial-model-obj)
$(CC) $(LFLAGS) -o $@ $(sd-radial-model-obj)
sd-opacity-src = \
$(SDLIB)/fluxmodel.c \
$(SIM5LIB)/sim5kerr.c \
sd-opacity.c
sd-opacity-obj = $(sd-opacity-src:.c=.o)
sd-opacity: $(sd-opacity-obj)
$(CC) $(LFLAGS) -o $@ $(sd-opacity-obj)
//************************************************************************
// sd-opacity.c -
//------------------------------------------------------------------------
// Written by:
// Michal Bursa (bursa@astro.cas.cz)
// Astronomical Institute
// Bocni II 1401/1a, 141-31 Praha 4, Czech Republic
//------------------------------------------------------------------------
// Created: 10.12.2012
//************************************************************************
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "fluxmodel.h"
#include "flux-sd.h"
#include "flux-nt.h"
#include "sim5kerr.h"
#include "sim5const.h"
#include "sim5math.h"
int main(int argc, char *argv[])
//***************************************************
{
if (argc!=6) {
fprintf(stderr, "Use: %s spin L alpha rmax nrad\n", argv[0]);
exit(-1);
}
double M = 10.0; // BH mass
double a = atof(argv[1]); // BH spin
double lum = atof(argv[2]); // BH luminosity
double alpha = atof(argv[3]); // alpha parameter
double rmax = atof(argv[4]); // maximal radius to consider
double rmin = r_horizon(a)*1.01; // minimal radius
int nrad = atoi(argv[5]); // number of rings
// initialize flux model
flux_model_init("./flux-sd.so", M, a, lum, alpha, NULL, SD_INIT_OPT_LUMINOSITY);
int ir;
for (ir=0; ir<nrad; ir++) {
double R = exp(log(rmin) + (log(rmax)-log(rmin))*(double)(ir)/(double)(nrad-1));
double T = flux_model_eval_temp(R); // [K]
double S = flux_model_eval_sigma(R); // [g cm^-2]
double H = flux_model_eval_h(R)*M*grav_radius_cgs; // [cm]
double l = flux_model_eval_l(R); // [?]
double V = flux_model_eval_vr(R); // [(c)]
double C = flux_model_eval_custom(R, SD_EVAL_TEMP_CENTRAL); // [K]
// central density
// density profile: // rho = rho0 * (1 - z^2/H^2)^N; N=3...polytropic index
// Sigma = \int_0^H rho dz = 16/35*H*rho0 (for N=3) ---> rho0 = 35/16*Sigma/H
double D = 35./16.*S/H; // [g cm^-3]
double tau_ff = 6.4e22 * 5.*M_PI/32. * pow(D,2.) * pow(C,-3.5) * H;
double tau_es = 0.34 * 16./35. * D * H;
double tau_eff = sqrt(tau_ff*(tau_ff+tau_es));
double tau_tot = tau_ff+tau_es;
double t_thermal = (H/tau_eff)/speed_of_light_cgs;
double t_freefall = (sqrt(2.*sqr3(R))/3. - 4./3.)*M*grav_radius_cgs/speed_of_light_cgs;
fprintf(stdout, "%.4e %e %e %e %e %e %e %e %e %e\n", R, C, T, D, H, tau_tot, tau_eff, tau_ff, tau_es, t_thermal/t_freefall);
}
flux_model_done();
return 0;
}
//************************************************************************
// sd-radial-model.c - print out radial model of disk
//------------------------------------------------------------------------
// Written by:
// Michal Bursa (bursa@astro.cas.cz)
// Astronomical Institute
// Bocni II 1401/1a, 141-31 Praha 4, Czech Republic
//------------------------------------------------------------------------
// Created: 18.10.2012
//************************************************************************
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "fluxmodel.h"
#include "flux-sd.h"
#include "flux-nt.h"
#include "sim5kerr.h"
#include "sim5const.h"
#include "sim5math.h"
typedef struct {
double R;
double T;
double S;
double H;
double l;
double V;
double Q;
double totL;
} radial_structure;
radial_structure* RS = NULL;
int nRS = 0;
int main(int argc, char *argv[])
//***************************************************
{
if (argc!=6) {
fprintf(stderr, "Use: %s spin L alpha rmax nrad\n", argv[0]);
exit(-1);
}
double M = 10.0; // BH mass
double a = atof(argv[1]); // BH spin
double lum = atof(argv[2]); // BH luminosity
double alpha = atof(argv[3]); // alpha parameter
double rmax = atof(argv[4]); // maximal radius to consider
double rmin = r_horizon(a)*1.05; // minimal radius
int nrad = atoi(argv[5]); // number of rings
// initialize flux model
flux_model_init("./flux-sd.so", M, a, lum, alpha, NULL, SD_INIT_OPT_LUMINOSITY);
nRS = nrad*100;
RS = (radial_structure*) calloc(nRS, sizeof(radial_structure));
double R0 = r_horizon(a);
// evaluate local model parameters at r=r0
int ir;
for (ir=0; ir<nRS; ir++) {
RS[ir].R = exp(log(rmin) + (log(rmax)-log(rmin))*(double)(ir)/(double)(nRS-1));
RS[ir].T = flux_model_eval_temp(RS[ir].R);
RS[ir].S = flux_model_eval_sigma(RS[ir].R);
RS[ir].H = flux_model_eval_h(RS[ir].R);
RS[ir].l = flux_model_eval_l(RS[ir].R);
RS[ir].V = flux_model_eval_vr(RS[ir].R);
RS[ir].totL = (ir>0?RS[ir-1].totL:0.0) + pow(RS[ir].T,4.)*sqr3(RS[ir].R)*(RS[ir].R-R0); // dL = F*r*dr
// vertical gravity (see Zhu et al. 2012)
double g00, g11, g22, g33, g03;
kerr_metric(a, RS[ir].R, 0.0, &g00, &g11, &g22, &g33, &g03);
double Omega = Omega_from_ell(RS[ir].l, g00, g33, g03);
double u[4];
u[0] = sqrt(-1.0/(g00 + 2.*Omega*g03 + sqr(Omega)*g33));
u[1] = 0.0;
u[2] = 0.0;
u[3] = u[0]*Omega;
if (fabs(RS[ir].V)>0.0) {
// radial velocity (V) is measured in fluid co-rotating frame
double e0[4], e1[4], e2[4], e3[4];
ortho_tetrad_U_azimuthal_motion(u, g00, g11, g22, g33, g03, e0, e1, e2, e3);
double u_loc[4];
u_loc[0] = sqrt(1.0+sqr(RS[ir].V)); // normalization to u.u=-1
u_loc[1] = RS[ir].V; // r-vector of the ortonorm frame points outwards & vr is negative
u_loc[2] = 0.0;
u_loc[3] = 0.0;
on2bl(u_loc, u, g00, g11, g22, g33, g03, e0, e1, e2, e3);
}
double u_t = u[0]*g00 + u[3]*g03;
double u_f = u[0]*g03 + u[3]*g33;
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;
R0 = RS[ir].R;
}
double dL = (log(RS[nRS-1].totL) - log(RS[0].totL))/(nrad-1);
int ID = 0;
for (ir=0; ir<nRS; ir++) {
if ((ir==nRS-1) || ((log(RS[ir].totL)-log(RS[0].totL)) > ID*dL)) {
fprintf(
stdout, "%04d %.4e %e %e %e %e %e %e\n",
ID+1, RS[ir].R, RS[ir].T, RS[ir].S, RS[ir].Q, RS[ir].H, RS[ir].l, RS[ir].V
);
ID++;
}
}
free(RS);
// close flux model
flux_model_done();
return 0;
}
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