/* * KING.C - tabulate a King (1966) Model for use with mkommod * * mar-91 Created Patrick Frisch (AdM) * jun-97 NEMO2 updates, more scaling options PJT/JJM * * Theory : (to be continued ...) * */ #include /* std NEMO include */ #include #include string defv[] = { "out=???\n Output File (required)", "w0=1\n Dimensionless central potential", "mass=1\n Total Mass of model", "rc=1\n Core Radius of model", "nsteps=101\n Number of steps to take", "epsilon=1e-6\n Accuracy for Integration", "tab=\n tabulate in an ASCII - file", "g=1\n Gravitational Constant", "VERSION=V1.0\n 27-jun-97 PJT", NULL, }; string usage = "tabulate a King (1966) Model for use with mkommod"; /* * Global definitions */ #define TABMAX 1000 #define TRYSTEPS 20 #define MAXEQS 10 /* * Global Parameters for the integrator */ local int nsteps, nequations, nok, nbad ; local real stepsize, stepmin, epsilon; int kmax, kount; /* in case a history record of */ double dxsav, xp[100], yp[10][100]; /* integration is needed */ /* used in NumRec functions */ /* * Global Parameters for King models */ local real radius[TABMAX],density[TABMAX],intmass[TABMAX], potential[TABMAX], distribution[TABMAX]; local real w0, rho0; /* * Function definitions ( I love prototyping (!!!)) */ int (*odefunc)(); /* This one will be exported to odeint */ local int function(real x, real *y, real *dxdy); local void initode(void); local void king(void ); local void odestep(real result[], real start, real end); local void out(void ); local real rho(real w); extern int rkqc(double*,double*,int,double*,double,double,double*, double*,double*,iproc); extern int odeint(double *,int,double,double,double,double,double, int*,int*,iproc,iproc); void nemo_main() { king(); out(); } local void king() { int i; real startpoint, endpoint, wstep, result[MAXEQS]; real sigma, rhoc, rc, mass, mscale, distr0, g; g = getdparam("g"); /* gravitational constant */ initode(); /* Potential Section : */ w0 = getdparam("w0"); /* First compute the unscaled model (depending on w0 only */ w0 = -ABS(w0); /* Be sure having negative potential */ rho0 = rho(w0); /* density unit */ potential[0] = w0; wstep = w0 / (nsteps-1); for(i = 1; i < nsteps; i++) potential[i] = potential[i-1] - wstep; /* Initial Value Section : */ radius[0] = 0.0; /* Starting Integration at center */ density[0] = 1.0; /* Starting density is central density */ intmass[0] = 0.0; result[0] = 0.0; /* Initialize vector for odeint with starting values */ result[1] = 2.0 / 3.0; /* need some finite value */ startpoint = w0; stepsize = (potential[1]-startpoint) / TRYSTEPS; /* give an initial estimate of integration stepsize */ stepmin = 0.0; odefunc = function; /* Stepping Section : */ for(i = 1; i < nsteps; i++) { endpoint = potential[i]; /* isn't necessary to define extra variables for start- and endpoint but makes it more readable */ odestep(result,startpoint,endpoint); /* main worker routine, calls odeint (see "Numerical Recipes") */ radius[i] = sqrt(result[0]); density[i] = rho(endpoint) / rho0; intmass[i] = 2. * result[0]*sqrt(result[0]) / result[1]; startpoint = endpoint; } printf("\n Unscaled model created !\n"); /* Scaling Section : */ /* * Calculate scaling factors from given total mass and core radius * * (see Binney&Tremaine [1987] for details) * */ mass = getdparam("mass"); /* total mass */ rc = getdparam("rc"); /* core radius */ mscale = mass / intmass[nsteps-1]; sigma = sqrt(g * mscale / rc); rhoc = 9.0 * sqr(sigma) / (FOUR_PI * g * sqr(rc)); distr0 = 9.0 / (2.0*FOUR_PI * g * sqr(rc) * sigma * sqrt(2.0) * rho0); /* This is to check the computed values against previously tabulated values */ printf("\nKing Model w0=%f M,rc,sigma,rhoc= %f %f %f %f\n", w0, mass, rc, sigma, rhoc); /* Tabulate Section : */ /* * Create final tables * */ for(i = 0;i < nsteps; i++) { radius[i] *= rc; density[i] *= rhoc; intmass[i] *= sqr(sigma) * rc / g ; /* NOTE : calculation of theoretical isotropic DF */ /* --------- */ /* for creation of anisotropic models see anisot.c */ distribution[i] = distr0 * (exp(-potential[i]) - 1); potential[i] *= sqr(sigma); dprintf(1,"\nRadius,Density,Mass,Pot,Dist= %f %f %f %f %f",radius[i], density[i],intmass[i],potential[i],distribution[i]); } } local void initode() { nequations = 2; /* the second order DE will be transformed to two first order DE's */ odefunc = function; /* Setup the function pointer */ epsilon = getdparam("epsilon"); /* for accuracy reasons */ nsteps = getiparam("nsteps"); } local void odestep(real y[], real start, real end) { int i; real ywork[MAXEQS]; for(i = 0; i < nequations; i++) ywork[i] = y[i]; odeint(ywork, nequations, start, end, epsilon, stepsize, stepmin, &nok, &nbad, function, rkqc); for(i = 0; i < nequations; i++) y[i] = ywork[i]; } local void out() { real ra=1; int i; string tabfile ; stream tab, outstr; if ( !streq( (tabfile = getparam("tab")), "" ) ) { tab = stropen(tabfile,"w"); fprintf(tab,"# Radius,Density,Mass,Potential,Distribution =\n"); for(i = 0; i < nsteps; i++) fprintf(tab,"%g %g %g %g %g\n",radius[i],density[i],intmass[i], potential[i],distribution[i]); } outstr = stropen(getparam("out"),"w"); /* This code comes directly from plummer.c */ put_set(outstr, "OsipkovMerrittModel"); put_data(outstr, "AnisoRadius", RealType, &ra, 0); put_data(outstr, "Ntab", IntType, &nsteps, 0); put_data(outstr, "Radius", RealType, radius, nsteps, 0); put_data(outstr, "Density", RealType, density, nsteps, 0); put_data(outstr, "Mass", RealType, intmass, nsteps, 0); put_data(outstr, "Potential", RealType, potential, nsteps, 0); put_data(outstr, "DistFunc", RealType, distribution, nsteps, 0); put_tes(outstr, "OsipkovMerrittModel"); strclose(outstr); } /* * The DE to be solved is (see King, A.J. 71, 64 [1966]): * * -x x" + 3/2 (x')^2 - 9/4 (rho(w)/rho0) (x')^3 = 0 * * where x = radius^2 and x' = dx/dw * * x, x' are stored in the two components of the vector y * For the numerical solution the DE is split in two first-order DE's: * * 1. dx/dw = z; * * 2. dz/dw = 3/2 z^2 (1 - 3/2 rho(w)/rho0 z) / x * * where z = z(w) = x' = y[1]; * */ local int function(real w, real y[], real dydw[]) { real z; dydw[0] = z = y[1]; if(y[0] > 1e-8) dydw[1] = 1.5*z*z*(1.0-1.5*rho(w)/rho0*z) / y[0]; else /* use an approximation near x=0 */ dydw[1] = 0.4*(1.0-pow(-w0,1.5)/rho0); return 0; } /* * The density in terms of w is given by: (see B&T, pp. 232) * * rho(w) = 0.5 * sqrt(PI) * exp(-w) * erf(y) - y - 2 * y^3 / 3 * * where y = sqrt( -w ) * Scaling factors will be neglected and handled later on */ local real rho(real w) { real y; real coeff=0.8862269254527579; /* sqrt(PI) / 2 */ if(w > 0.0) return 0.0; y = sqrt(-w); return exp(-w)*coeff*erf(y) - y - 2.0*pow(y,3.0)/3.0; }