-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathballic.h
140 lines (119 loc) · 4.76 KB
/
ballic.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
/*
** The header file for ballic.
*/
#ifndef BALLIC_HINCLUDED
#define BALLIC_HINCLUDED
#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;
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;
int nLayer;
/*
** This array contains the material number for each layer (e.g., IRON, GRANITE).
*/
int *iLayer;
double *MLayer; // Mass of each layer
/*
** 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;
int *mat;
double dr;
double R;
} MODEL;
const double fact = 1.0;
/* Icosahedron. */
ICOSA *icosaInit(void);
void icosaPix2Vec(ICOSA *ctx,int i,int resolution,double *vec);
/* HEALPix. */
void pix2vec_ring( long nside, long ipix, double *vec);
MODEL *modelInit(double M,double ucore);
double drhodr(MODEL *model,int iLayer,double r,double rho,double M,double u);
double dudr(MODEL *model,int iLayer,double r,double rho,double M,double u);
double dudrho(MODEL *model,int iLayer,double rho,double u);
double dMdr(double r,double rho);
void modelSolveBC(MODEL *model, double *prho, double *pu, int iLayer1, int iLayer2);
void modelSolveComponent(MODEL *model,int iLayer,int bSetModel,int bLastLayer,int *pIndex,double h,double *prho1,double *pu1,double *pM1,double M2,double *pR);
double modelSolveTwoComponent(MODEL *model,int bSetModel,double rho,double u,double h,double *pR);
void modelWriteToFile(MODEL *model);
double modelSolveAll(MODEL *model,int bSetModel,double rhoc,double uc,double h,double *pR);
double modelSolve(MODEL *model,double M);
double MLookup(MODEL *model,double r);
double rhoLookup(MODEL *model,double r);
double uLookup(MODEL *model,double r);
/* Basic ballic routines. */
double Fzero(MODEL *model,int bIcosa,double r,double ri,double m,int ns);
double rShell(MODEL *model,double m,double ri);
double rShell2(MODEL *model,int bIcosa,double m,double ri,int ns);
#endif