https://github.com/chreinhardt/ballic
Revision 5dc1ea022a08d9cc896e5dca450ab5b95036c364 authored by Christian Reinhardt on 02 January 2023, 10:04:03 UTC, committed by Christian Reinhardt on 02 January 2023, 10:04:03 UTC
1 parent 13d71d8
Tip revision: 5dc1ea022a08d9cc896e5dca450ab5b95036c364 authored by Christian Reinhardt on 02 January 2023, 10:04:03 UTC
Update code for the new format of the model files.
Update code for the new format of the model files.
Tip revision: 5dc1ea0
ballic.single.c
/*
* Copyright (c) 2014-2019 Joachim Stadel and Christian Reinhardt.
*
* ballic provides a low noise particle representation of equilibrium models.
*
* The single version builds a 1D equilibrium model for a given material and boundary conditions and
* then creates the particle representation which is stored in ballic.std.
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include "tipsydefs/tipsy.h"
#include "tillotson/tillotson.h"
#define max(A,B) ((A) > (B) ? (A) : (B))
#define min(A,B) ((A) > (B) ? (B) : (A))
// Isentropic thermal profile
#define BALLIC_U_ISENTROPIC
typedef struct icosa_struct {
float R[180];
float v[36];
} ICOSA;
ICOSA *icosaInit(void) {
ICOSA *ctx;
ctx = malloc(sizeof(ICOSA));
assert(ctx != NULL);
compute_matrices_(&ctx->R);
compute_corners_(&ctx->v);
return(ctx);
}
void icosaPix2Vec(ICOSA *ctx,int i,int resolution,double *vec) {
float v[3];
pixel2vector_(&i,&resolution,&ctx->R,&ctx->v,v);
vec[0] = v[0];
vec[1] = v[1];
vec[2] = v[2];
}
/* -----------------------------------------------------------------------------
*
* Copyright (C) 1997-2005 Krzysztof M. Gorski, Eric Hivon,
* Benjamin D. Wandelt, Anthony J. Banday,
* Matthias Bartelmann,
* Reza Ansari & Kenneth M. Ganga
*
*
* This file is part of HEALPix.
*
* HEALPix is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* HEALPix is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with HEALPix; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*
* For more information about HEALPix see http://healpix.jpl.nasa.gov
*
*----------------------------------------------------------------------------- */
void pix2vec_ring( long nside, long ipix, double *vec) {
/*
c=======================================================================
c gives theta and phi corresponding to pixel ipix (RING)
c for a parameter nside
c=======================================================================
*/
int nl2, nl4, npix, ncap, iring, iphi, ip, ipix1;
double fact1, fact2, fodd, hip, fihip;
double PI=M_PI;
double z, sz, phi;
// PARAMETER (pi = 3.1415926535897932384626434d0)
// parameter (ns_max = 8192) ! 2^13 : largest nside available
int ns_max=8192;
if( nside<1 || nside>ns_max ) {
fprintf(stderr, "%s (%d): nside out of range: %ld\n", __FILE__, __LINE__, nside);
exit(0);
}
npix = 12*nside*nside; // ! total number of points
if( ipix<0 || ipix>npix-1 ) {
fprintf(stderr, "%s (%d): ipix out of range: %ld\n", __FILE__, __LINE__, ipix);
exit(0);
}
ipix1 = ipix + 1; // in {1, npix}
nl2 = 2*nside;
nl4 = 4*nside;
ncap = 2*nside*(nside-1);// ! points in each polar cap, =0 for nside =1
fact1 = 1.5*nside;
fact2 = 3.0*nside*nside;
if( ipix1 <= ncap ) { //! North Polar cap -------------
hip = ipix1/2.;
fihip = floor(hip);
iring = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;// ! counted from North pole
iphi = ipix1 - 2*iring*(iring - 1);
z = 1. - iring*iring / fact2 ;
phi = (1.*iphi - 0.5) * PI/(2.*iring);
}
else if( ipix1 <= nl2*(5*nside+1) ) {//then ! Equatorial region ------
ip = ipix1 - ncap - 1;
iring = (int)floor( ip / nl4 ) + nside;// ! counted from North pole
iphi = (int)fmod(ip,nl4) + 1;
fodd = 0.5 * (1 + fmod((double)(iring+nside),2));// ! 1 if iring+nside is odd, 1/2 otherwise
z = (nl2 - iring) / fact1;
phi = (1.*iphi - fodd) * PI /(2.*nside);
}
else {//! South Polar cap -----------------------------------
ip = npix - ipix1 + 1;
hip = ip/2.;
/* bug corrige floor instead of 1.* */
fihip = floor(hip);
iring = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;// ! counted from South pole
iphi = (int)(4.*iring + 1 - (ip - 2.*iring*(iring-1)));
z = -1. + iring*iring / fact2 ;
phi = (1.*iphi - 0.5) * PI/(2.*iring);
}
sz = sqrt( 1.0 - z*z );
vec[0] = sz * cos(phi);
vec[1] = sz * sin(phi);
vec[2] = z;
}
double Packed49[49][3] = {
{1.0820379E-008,-3.2300459E-006, 0.7790824},
{ 0.0851419, -0.0604194, 0.3497294},
{-0.0914288, 0.4137149, 0.6537967},
{-0.4233704, -0.0167043, 0.6537949},
{ 0.4161716, 0.0795128, 0.6537953},
{ 0.2606476, -0.3340458, 0.6537932},
{-0.1783143, -0.3843534, 0.6537934},
{ 0.3214161, 0.4881783, 0.5151146},
{-0.4918671, 0.3793286, 0.4702615},
{-0.0929470, 0.3254455, 0.2208710},
{-0.5385916, -0.3977282, 0.3983726},
{ 0.6544342, -0.1767253, 0.3839966},
{-0.0464433, -0.6884808, 0.3616719},
{ 0.6543453 , 0.2637903 , 0.3304789},
{ 0.3839234, -0.5971428, 0.3209247},
{-0.7108776, 0.0012718, 0.3187802},
{-0.1876019, -0.3352532, 0.1368439},
{ 0.0704146, 0.7351211, 0.2481284},
{-0.3638534, 0.6695231, 0.1622308},
{ 0.3033485, 0.2022102, 0.0692749},
{ 0.4742969, 0.6067250, 0.1178835},
{-0.6795831, 0.3744220, 0.0703152},
{-0.3794767, 0.0482692, 0.0303653},
{ 0.6538247, -0.4232979, 0.0173632},
{ 0.2332925, -0.2885923, 0.0015460},
{ 0.7790813, -5.7994478E-013, -0.0012951},
{-0.5429419, -0.5585797, -0.0131201},
{-0.1452212, -0.7628716, -0.0625065},
{-0.7541588, -0.1768841, -0.0832219},
{ 0.2928920, -0.7156702, -0.0948672},
{-0.0815266, 0.7567586, -0.1662504},
{ 0.6534047, 0.3539813, -0.2339385},
{-0.4662671, 0.5642231, -0.2668646},
{ 0.3250845, 0.6429856, -0.2964101},
{-0.6822678, 0.1837015, -0.3282282},
{ 0.4930282, -0.4578927, -0.3927173},
{-0.3428409, -0.5690840, -0.4069064},
{ 0.6530941, -0.0471244, -0.4221572},
{-0.1036092, 0.2867325, -0.2191358},
{ 0.0858176, -0.5795982, -0.5134887},
{-0.5513357, -0.1947856, -0.5148367},
{ 0.0032634, 0.5253046, -0.5753380},
{-0.1660916, -0.1851457, -0.2781839},
{ 0.3945018, 0.3205518, -0.5904102},
{-0.3693146, 0.2921726, -0.6206539},
{ 0.2268184, 0.0121013, -0.3221586},
{ 0.2750635, -0.2203113, -0.6948182},
{-0.1633231, -0.2729136, -0.7112054},
{ 0.0133638, 0.1279994, -0.7683794}
};
typedef struct model_ctx {
/* Material coefficients from the Tillotson EOS. */
TILLMATERIAL *tillMat;
/*
** Some unit conversion factors.
*/
double dKpcUnit;
double dMsolUnit;
/*
** The lookup table for the equilibrium model.
*/
int nTableMax;
int nTable;
double uc; /* u at r = 0 */
double *M;
double *rho;
double *u;
double *r;
double dr;
double R;
} MODEL;
MODEL *modelInit(double ucore, int iMat) {
/* Initialize the model */
MODEL *model;
model = malloc(sizeof(MODEL));
assert(model != NULL);
model->dKpcUnit = 2.06701e-13;
model->dMsolUnit = 4.80438e-08;
model->tillMat = malloc(sizeof(TILLMATERIAL));
/* Check if the Tillotson library has the right version. */
if (TILL_VERSION_MAJOR != 3) {
fprintf(stderr, "modelInit: Tillotson library has the wrong version (%s)\n", TILL_VERSION_TEXT);
exit(1);
}
fprintf(stderr, "Tillotson EOS library version: %s\n", TILL_VERSION_TEXT);
fprintf(stderr, "\n");
/*
* Initialize one material:
*
* i=0: Ideal gas
* i=1: Granite
* i=2: Iron
* i=3: Basalt
* i=4: Ice
*
* All materials are defined in tillotson.h.
*/
model->tillMat = tillInitMaterial(iMat, model->dKpcUnit, model->dMsolUnit);
tillPrintMat(model->tillMat);
/* model->uFixed = uFixed/model->dErgPerGmUnit; */
model->uc = ucore;
model->nTableMax = 10000;
model->M = malloc(model->nTableMax*sizeof(double));
assert(model->M != NULL);
model->rho = malloc(model->nTableMax*sizeof(double));
assert(model->rho != NULL);
model->u = malloc(model->nTableMax*sizeof(double));
assert(model->u != NULL);
model->r = malloc(model->nTableMax*sizeof(double));
assert(model->r != NULL);
model->dr = 0.0;
model->nTable = 0;
return(model);
}
double drhodr(MODEL *model,double r,double rho,double M,double u);
double dudr(MODEL *model,double r,double rho,double M,double u);
/*
* Calculate drhodr to solve for the equilibrium model.
*/
double drhodr(MODEL *model,double r,double rho,double M,double u) {
double dPdrho,dPdu;
dPdrho=tilldPdrho(model->tillMat, rho, u); // dP/drho at u=const.
dPdu = tilldPdu(model->tillMat, rho, u);; // dP/du at rho=const.
/*
* drho/dr = -G*M*rho/(dPdrho+dPdu*dudrho)
*/
assert(r >= 0.0);
if (r > 0.0) {
// We assume G=1
return(-M*rho/(r*r*(dPdrho + dPdu*tilldudrho(model->tillMat,rho,u))));
}
else {
return(0.0);
}
}
/*
* We assume that the thermal profile is isentropic.
*/
double dudr(MODEL *model,double r,double rho,double M,double u) {
return(tilldudrho(model->tillMat, rho, u)*drhodr(model, r, rho, M, u));
}
/*
* This derivative is independent of the model and only involves geometry.
*/
double dMdr(double r,double rho) {
assert(r >= 0.0);
return(4.0*M_PI*r*r*rho);
}
const double fact = 1.0;
/*
* Calculate the gravitational accelleration due to the enclosed mass M(R).
*/
double CalcGrav(double r, double M)
{
assert(r>=0.0);
if (r > 0.0)
{
// Again assuming G=1
return(-M/(r*r));
} else {
return(0.0);
}
}
/*
* This function solves the model as an initial value problem with rho_initial = rho and
* M_initial = 0 at r = 0. This function returns the mass when rho == model->tillMat[i]->rho0.
*/
double midPtRK(MODEL *model,int bSetModel,double rho,double h,double *pR) {
FILE *fp;
double M = 0.0;
double r = 0.0;
double u = model->uc;
double k1rho,k1M,k1u,k2rho,k2M,k2u,x;
int i;
if (bSetModel) {
i = 0;
model->rho[i] = rho;
model->M[i] = M;
model->u[i] = u;
model->r[i] = r;
fp = fopen("ballic.model","w");
assert(fp != NULL);
fprintf(fp,"%g %g %g %g %g\n",r,rho,M,u,CalcGrav(r,M));
++i;
}
while (rho > fact*model->tillMat->rho0) {
/*
** Midpoint Runga-Kutta (2nd order).
*/
k1rho = h*drhodr(model,r,rho,M,u);
k1M = h*dMdr(r,rho);
k1u = h*dudr(model,r,rho,M,u);
k2rho = h*drhodr(model,r+0.5*h,rho+0.5*k1rho,M+0.5*k1M,u+0.5*k1u);
k2M = h*dMdr(r+0.5*h,rho+0.5*k1rho);
k2u = h*dudr(model,r+0.5*h,rho+0.5*k1rho,M+0.5*k1M,u+0.5*k1u);
rho += k2rho;
M += k2M;
u += k2u;
r += h;
if (bSetModel) {
model->rho[i] = rho;
model->M[i] = M;
model->u[i] = u;
model->r[i] = r;
fprintf(fp,"%g %g %g %g %g\n",r,rho,M,u,CalcGrav(r,M));
++i;
}
}
/*
** Now do a linear interpolation to rho == fact*rho0.
*/
x = (fact*model->tillMat->rho0 - rho)/k2rho;
assert(x <= 0.0);
r += h*x;
M += k2M*x;
rho += k2rho*x;
u += k2u*x;
if (bSetModel) {
--i;
model->M[i] = M;
model->r[i] = r;
model->rho[i] = rho;
model->u[i] = u;
fprintf(fp,"%g %g %g %g %g\n",r,rho,M,u,CalcGrav(r,M));
fclose(fp);
++i;
model->nTable = i;
model->dr = h;
}
*pR = r;
return(M);
}
/*
* Solve the 1D equilibrium model for the desired mass.
*/
double modelSolve(MODEL *model, double M) {
const int nStepsMax = 10000;
int bSetModel;
double rmax;
double dr,R;
double a,Ma,b,Mb,c,Mc;
/*
* First estimate the maximum possible radius.
*/
R = cbrt(3.0*M/(4.0*M_PI*model->tillMat->rho0));
dr = R/nStepsMax;
a = 1.01*model->tillMat->rho0; /* starts with 1% larger central density */
Ma = midPtRK(model,bSetModel=0,a,dr,&R);
fprintf(stderr,"first Ma:%g R:%g\n",Ma,R);
b = a;
Mb = 0.5*M;
while (Ma > M) {
b = a;
Mb = Ma;
a = 0.5*(model->tillMat->rho0 + a);
Ma = midPtRK(model,bSetModel=0,a,dr,&R);
}
while (Mb < M) {
b = 2.0*b;
Mb = midPtRK(model,bSetModel=0,b,dr,&R);
fprintf(stderr,"first Mb:%g R:%g\n",Mb,R);
}
// (CR) Debug
fprintf(stderr,"Root bracketed.\n");
/*
* Root bracketed by (a,b).
*/
while (Mb-Ma > 1e-10*Mc) {
c = 0.5*(a + b);
Mc = midPtRK(model,bSetModel=0,c,dr,&R);
if (Mc < M) {
a = c;
Ma = Mc;
}
else {
b = c;
Mb = Mc;
}
// fprintf(stderr,"c:%.10g Mc:%.10g R:%.10g\n",c/model->tillMat[0]->rho0,Mc,R);
}
/*
* Solve it once more setting up the lookup table.
*/
fprintf(stderr,"rho_core: %g cv: %g uc: %g (in system units)\n",c,model->tillMat->cv,model->uc);
Mc = midPtRK(model,bSetModel=1,c,dr,&R);
model->R = R;
return c;
}
double MLookup(MODEL *model,double r) {
double x,xi,dr;
int i;
i = model->nTable-1;
if (r >= model->r[i]) return(model->M[i]*(1.0 + log(r-model->r[i]+1)));
x = r/model->dr;
xi = floor(x);
assert(xi >= 0.0);
x -= xi;
i = (int)xi;
if (i < 0) {
fprintf(stderr,"ERROR r:%.14g x:%.14g xi:%.14g i:%d\n",r,x,xi,i);
}
assert(i >= 0);
if (i < model->nTable-2) {
return(model->M[i]*(1.0-x) + model->M[i+1]*x);
}
if (i == model->nTable-2) {
dr = model->r[i+1] - model->r[i];
x = r/dr;
xi = floor(x);
x -= xi;
return(model->M[i]*(1.0-x) + model->M[i+1]*x);
}
else {
i = model->nTable - 1;
return(model->M[i]*(1.0 + log(r-model->r[i]+1)));
}
}
double rhoLookup(MODEL *model,double r) {
double x,xi,dr;
int i;
i = model->nTable-1;
if (r >= model->r[i]) return(model->rho[i]*exp(-(r-model->r[i])));
x = r/model->dr;
xi = floor(x);
assert(xi >= 0.0);
x -= xi;
i = (int)xi;
if (i < 0) {
fprintf(stderr,"ERROR r:%.14g x:%.14g xi:%.14g i:%d\n",r,x,xi,i);
}
assert(i >= 0);
if (i < model->nTable-2) {
return(model->rho[i]*(1.0-x) + model->rho[i+1]*x);
}
if (i == model->nTable-2) {
dr = model->r[i+1] - model->r[i];
x = r/dr;
xi = floor(x);
x -= xi;
return(model->rho[i]*(1.0-x) + model->rho[i+1]*x);
}
else {
i = model->nTable - 1;
return(model->rho[i]*exp(-(r-model->r[i])));
}
}
double uLookup(MODEL *model,double r) {
double x,xi,dr;
int i;
i = model->nTable-1;
if (r >= model->r[i]) return(model->u[i]*exp(-(r-model->r[i])));
x = r/model->dr;
xi = floor(x);
assert(xi >= 0.0);
x -= xi;
i = (int)xi;
if (i < 0) {
fprintf(stderr,"ERROR r:%.14g x:%.14g xi:%.14g i:%d\n",r,x,xi,i);
}
assert(i >= 0);
if (i < model->nTable-2) {
return(model->u[i]*(1.0-x) + model->u[i+1]*x);
}
if (i == model->nTable-2) {
dr = model->r[i+1] - model->r[i];
x = r/dr;
xi = floor(x);
x -= xi;
return(model->u[i]*(1.0-x) + model->u[i+1]*x);
}
else {
i = model->nTable - 1;
return(model->u[i]*exp(-(r-model->r[i])));
}
}
double Fzero(MODEL *model,int bIcosa,double r,double ri,double m,int ns) {
long npix = (bIcosa)?(40*ns*(ns-1)+12):(12*ns*ns);
return(MLookup(model,r)-MLookup(model,ri)-npix*m);
}
double rShell(MODEL *model,double m,double ri) {
double a = ri;
double b = 1.0;
double c,Mc,Ma;
Ma = MLookup(model,a);
Mc = 0.0;
while (m > (MLookup(model,b)-Ma)) b *= 2.0;
while (fabs(m-(Mc-Ma)) > 1e-10*m) {
c = 0.5*(a + b);
Mc = MLookup(model,c);
// fprintf(stderr,"c:%.7g Mc:%.7g\n",c,Mc);
if (m > (Mc-Ma)) a = c;
else b = c;
}
return c;
}
double rShell2(MODEL *model,int bIcosa,double m,double ri,int ns) {
double a = ri;
double b = 1.0;
double c;
double z;
z = 1.0;
while (Fzero(model,bIcosa,b,ri,m,ns) < 0) b *= 2.0;
while ((b-a)/c > 1e-10) {
c = 0.5*(a + b);
z = Fzero(model,bIcosa,c,ri,m,ns);
if (z < 0) a = c;
else b = c;
// printf("c:%.14g M(c):%.14g\n",c,MLookup(model,c));
}
return c;
}
void main(int argc, char **argv) {
const int bCentral = 1;
int bIcosa = 0;
const int bRandomRotate = 1;
TCTX out;
long ns,npix,ipix,na,nb;
struct gas_particle gp;
double r[3];
double rhoCenter;
double ri,ro,rs,roa,rob,ros,rta,rtb,rts;
double m,mTot,l1,l2,nsf;
double x,y,ang1,ang2,ang3;
int j,iShell,nDesired,nReached,nLast;
double theta,phi;
float rr[3];
ICOSA *ctx;
int nShell,nMaxShell;
double *rsShell;
long *nsShell;
long *isShell;
double *xyz;
// Model
MODEL *model;
double ucore;
int iMat;
int iter;
int nSmooth;
double d,rho,u,eta,w0,dPdrho;
FILE *fpi,*fpo;
int iRet,i;
/*
** These variables are used to find the optimal softening.
*/
double l1max = 0.0;
double l2max = 0.0;
if (argc != 5) {
fprintf(stderr,"Usage: ballic <nDesired> <TotalMass> <ucore> <iMat> >myball.std\n");
exit(1);
}
nDesired = atoi(argv[1]);
mTot = atof(argv[2]);
ucore = atof(argv[3]);
iMat = atoi(argv[4]);
/* Assure that the parameters make sense. */
assert(nDesired > 0);
assert(mTot > 0.0);
assert(ucore > 0.0);
assert(iMat >= 0);
model = modelInit(ucore, iMat);
rhoCenter = modelSolve(model,mTot);
m = mTot/nDesired; /* a first guess at the particle mass */
/*
** Initialize icosahedral parameters.
*/
ctx = icosaInit();
#if (0)
/*
** Test the icosahedral grid.
*/
ns = 4;
npix = 40*ns*(ns-1) + 12;
for (ipix=0;ipix<nipix;++ipix) {
icosaPix2Vec(ctx,ipix,ns,r);
printf("%d %g %g %g\n",i,r[0],r[1],r[2]);
}
#endif
/*
** Initialize the array of shell radii and resolutions.
*/
nMaxShell = 1000;
nShell = 0;
rsShell = malloc(nMaxShell*sizeof(double));
assert(rsShell != NULL);
nsShell = malloc(nMaxShell*sizeof(long));
assert(nsShell != NULL);
isShell = malloc(nMaxShell*sizeof(long));
assert(nsShell != NULL);
/*
** Phase 1: We first determine the number of particles per shell such that
** the ratio of radial to tangential lengths of their volumes is as close
** to 1 as possible. This fixes the total number of particles in the model.
*/
fprintf(stderr,"PHASE 1: ODE approach\n");
for (iter=0;iter<2;++iter) {
nReached = 0;
if (bCentral) {
ro = rShell(model,m,0.0);
}
else {
ro = 0.0;
}
iShell = 0;
for (;;) {
ri = ro;
na = 1;
roa = rShell2(model,bIcosa,m,ri,na);
npix = (bIcosa)?(40*na*(na-1)+12):(12*na*na);
l1 = roa-ri;
l2 = sqrt(M_PI/npix)*(roa+ri);
rta = l1/l2;
nb = 16;
do {
nb *= 2;
rob = rShell2(model,bIcosa,m,ri,nb);
npix = (bIcosa)?(40*nb*(nb-1)+12):(12*nb*nb);
l1 = rob-ri;
l2 = sqrt(M_PI/npix)*(rob+ri);
rtb = l1/l2;
} while (rtb < 1.0);
while (nb - na > 1) {
ns = (na+nb)/2;
ros = rShell2(model,bIcosa,m,ri,ns);
npix = (bIcosa)?(40*ns*(ns-1)+12):(12*ns*ns);
l1 = ros-ri;
l2 = sqrt(M_PI/npix)*(ros+ri);
rts = l1/l2;
/* fprintf(stderr,"ns:%d rts:%g\n",ns,rts); */
if (rts < 1.0) {
na = ns;
roa = ros;
rta = rts;
}
else {
nb = ns;
rob = ros;
rtb = rts;
}
}
/*
if (iShell >= 5) {
fprintf(stderr,"na:%d rta:%g\n",na,rta);
fprintf(stderr,"nb:%d rtb:%g\n",nb,rtb);
}
*/
/*
** if the two possible ratios differ by less that 1% then we favour
** the higher resolution spherical grid (nb).
*/
if (1/rta+0.01 < rtb) {
ro = roa;
ns = na;
rts = rta;
}
else {
ro = rob;
ns = nb;
rts = rtb;
}
npix = (bIcosa)?(40*ns*(ns-1)+12):(12*ns*ns);
if (iShell == nMaxShell) {
nMaxShell *= 2;
rsShell = realloc(rsShell,nMaxShell*sizeof(double));
assert(rsShell != NULL);
nsShell = realloc(nsShell,nMaxShell*sizeof(long));
assert(nsShell != NULL);
}
nsShell[iShell] = ns;
/* fprintf(stderr,"nReached:%d npix:%d\n",nReached,npix);*/
if ((nReached + npix) < nDesired) {
nReached += npix;
fprintf(stderr,"iShell:%d ns:%d radial/tangential:%g\n",iShell,ns,rts);
++iShell;
}
else {
nShell = iShell;
break;
}
} /* end of iShell loop */
ns = nsShell[iShell-1];
npix = (bIcosa)?(40*ns*(ns-1)+12):(12*ns*ns);
if (nDesired - nReached > npix/2) {
nReached += npix;
nsShell[nShell] = ns;
fprintf(stderr,"iShell:%d ns:%d radial/tangential:??? (added)\n",nShell,ns);
nShell++;
}
fprintf(stderr,"nReached:%d old mass:%.7g new mass:%.7g\n",
nReached,m,mTot/nReached);
m = mTot/nReached;
nDesired = nReached+1;
}
/*
** Phase 2: With the numbers of particles in each of the shells now fixed
** we continue by recalculating the radii of the shells based on the updated
** particle mass.
*/
fprintf(stderr,"PHASE 2\n");
if (bCentral) {
ro = rShell(model,m,0.0);
}
else {
ro = 0.0;
}
for (iShell=0;iShell<nShell;++iShell) {
ri = ro;
ns = nsShell[iShell];
ro = rShell2(model,bIcosa,m,ri,ns);
npix = (bIcosa)?(40*ns*(ns-1)+12):(12*ns*ns);
l1 = ro-ri;
l2 = sqrt(M_PI/npix)*(ro+ri);
rts = l1/l2;
/*
** Here we calculate the optimal softening for gravity using the
** largest bin size. l1: radial size, l2: tangental size
*/
if (l1 > l1max) {
l1max = l1;
}
if (l2 > l2max) {
l2max = l2;
}
if (bIcosa) {
d = 1.0;
if (iShell == 0) d = 2.5;
if (iShell == 1) d = 1.0;
if (iShell == 2) d = 3.0;
}
else {
d = 1.9;
if (iShell == 0) d = 2.0;
if (iShell == 1) d = 1.5;
}
rs = rsShell[iShell] = pow(0.5*(pow(ri,d) + pow(ro,d)),1/d);
rho = rhoLookup(model,rs);
u = uLookup(model,rs); /* We also have to look up u from a table */
eta = rho/model->tillMat->rho0;
/* This was the old code using a constant internal energy uFixed.
w0 = model->uFixed/(model->par.u0*eta*eta) + 1.0;
dPdrho = (model->par.a + (model->par.b/w0)*(3 - 2/w0))*model->uFixed +
(model->par.A + 2*model->par.B*(eta - 1))/model->par.rho0;
fprintf(stderr,"iShell:%d r:%g M:%g rho:%g ns:%d radial/tangential:%g dr:%g <? Jeans:%g Gamma:%g\n",iShell,rs,MLookup(model,rs),rho,ns,rts,ro-ri,sqrt(dPdrho/rho),Gamma(model,rho,model->uFixed));
*/
w0 = u/(model->tillMat->u0*eta*eta) + 1.0;
dPdrho = (model->tillMat->a + (model->tillMat->b/w0)*(3 - 2/w0))*u +
(model->tillMat->A + 2*model->tillMat->B*(eta - 1))/model->tillMat->rho0;
// fprintf(stderr,"iShell:%d r:%g M:%g rho:%g u:%g ns:%d radial/tangential:%g dr:%g <? Jeans:%g Gamma:%g\n",iShell,rs,MLookup(model,rs),rho,u,ns,rts,ro-ri,sqrt(dPdrho/rho),Gamma(model,rho,u));
}
/*
** Now generate the coordinates of all the particles in each shell as they are on the unit
** sphere. This simplifies the later adjusting of the radii of the shells (unit sphere
** coordinates are independent of this.
*/
/*
** Phase 3: With the masses and numbers of the particles fixed and the
** as well as their angular positions, we now attempt to use an SPH density
** estimate to calculate the mean radial pressure force on the shell
** and converge on the radius of the shell such that it is balanced by
** gravity on the shell (analytic).
*/
if (bCentral) {
/*
** First adjust the density of the central particle.
*/
}
/*
** Now output the particles.
*/
TipsyInitialize(&out,0,NULL);
for (j=0;j<3;++j) gp.vel[j] = 0.0;
gp.phi = 0.0;
/*
** As a first guess we use max(l1max,l2max) for the softening.
*/
fprintf(stderr,"l1max: %g l2max: %g epsilon: %g\n",l1max,l2max,MAX(l1max,l2max));
gp.hsmooth = MAX(l1max,l2max); /* is actually eps for gasoline */
gp.mass = m;
fprintf(stderr,"hsmooth=%g\n",gp.hsmooth);
gp.hsmooth=0.001;
//gp.temp = model->uFixed; /* Christian's version of gasoline uses thermal energy instead of temperature as input! */
nLast = nReached;
nReached = 0;
if (bCentral) {
for (j=0;j<3;++j) gp.pos[j] = 0.0;
gp.temp = uLookup(model, 0);
// Dont forget to set the material for the central particle
gp.metals = iMat;
TipsyAddGas(out,&gp);
}
for (iShell=0;iShell<nShell;++iShell) {
rs = rsShell[iShell];
ns = nsShell[iShell];
npix = (bIcosa)?(40*ns*(ns-1)+12):(12*ns*ns);
nReached += npix;
ang1 = 2.0*M_PI*rand()/(RAND_MAX+1.0);
ang2 = 2.0*M_PI*rand()/(RAND_MAX+1.0);
ang3 = 2.0*M_PI*rand()/(RAND_MAX+1.0);
/* Experiment with grav. softening. */
if (iShell < nShell-1)
{
// Inner shell
gp.hsmooth = 0.001;
} else {
// Most outer shell
gp.hsmooth = 0.01;
}
for (ipix = 0;ipix < npix;++ipix) {
if (bIcosa) icosaPix2Vec(ctx,ipix,ns,r);
else pix2vec_ring(ns,ipix,r);
if (bRandomRotate) {
y = r[1]*cos(ang1) - r[2]*sin(ang1);
r[2] = r[1]*sin(ang1) + r[2]*cos(ang1);
x = r[0]*cos(ang2) - r[2]*sin(ang2);
r[2] = r[0]*sin(ang2) + r[2]*cos(ang2);
r[0] = x;
r[1] = y*cos(ang3) - r[2]*sin(ang3);
r[2] = y*sin(ang3) + r[2]*cos(ang3);
}
for (j=0;j<3;++j) gp.pos[j] = rs*r[j];
// rho = rhoLookup(model,rs);
gp.temp = uLookup(model,rs);
// Save Material
gp.metals = iMat;
TipsyAddGas(out,&gp);
}
}
fprintf(stderr,"Writing %d particles. Model R:%g Last Shell r:%g\n",nReached,model->R,rsShell[nShell-1]);
TipsyWriteAll(out,0.0,"ballic.std");
TipsyFinish(out);
#if 0
/* Smooth with gather doesn work with my version! */
system("sleep 1; smooth -s 128 density <ballic.std; sleep 1");
fpi = fopen("smooth.den","r");
assert(fpi != NULL);
fpo = fopen("ballic.den","w");
assert(fpo != NULL);
iRet = fscanf(fpi,"%d",&nReached);
assert(iRet == 1);
if (bCentral) {
iRet = fscanf(fpi,"%lf",&rho);
assert(iRet == 1);
fprintf(fpo,"0.0 %g\n",rho);
nReached -= 1;
}
for (iShell=0;iShell<nShell;++iShell) {
rs = rsShell[iShell];
ns = nsShell[iShell];
npix = (bIcosa)?(40*ns*(ns-1)+12):(12*ns*ns);
nReached -= npix;
for (ipix = 0;ipix < npix;++ipix) {
iRet = fscanf(fpi,"%lf",&rho);
assert(iRet == 1);
fprintf(fpo,"%g %g\n",rs,rho);
}
}
fclose(fpi);
fclose(fpo);
assert(nReached == 0);
#endif
}
Computing file changes ...