Skip to content

Commit 849a02a

Browse files
ctmblctmbl
committed
Change eos_* signature to pass the whole star composition
Changing eos_freeeos signature forces to change eos_calc signatur and also, for consistency purposes, the other eos: eos_ideal, eos_idealrad, eos_opal functions signature it's even a bigger mess because atm_onelayer and so atm_calc calls eos_calc, so there siganture had to be changed too a bunch of TODO had been added too and will have to be addressed
1 parent 01f4986 commit 849a02a

8 files changed

Lines changed: 53 additions & 54 deletions

File tree

‎src/include/physics.h‎

Lines changed: 13 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -34,12 +34,12 @@ class composition_map : public matrix_map {
3434

3535
int opa_calc(const matrix &X,double Z,const matrix &T,const matrix &rho,
3636
opa_struct &opa);
37-
int eos_calc(const matrix &X,double Z,const matrix &T,const matrix &p,
38-
matrix &rho,eos_struct &eos);
37+
int eos_calc(const composition_map &, const matrix &T, const matrix &p,
38+
matrix &rho, eos_struct &eos);
3939
int nuc_calc(const composition_map &, const matrix &T, const matrix &rho,
4040
nuc_struct &nuc);
41-
int atm_calc(const matrix &X,double Z,const matrix &g,const matrix &Teff,
42-
const char *eos_name,const char *opa_name,atm_struct &atm);
41+
int atm_calc(const composition_map &, const matrix &g, const matrix &Teff,
42+
const char *eos_name, const char *opa_name, atm_struct &atm);
4343

4444
int opa_opal(const matrix &X,double Z,const matrix &T,const matrix &rho,
4545
opa_struct &opa);
@@ -58,17 +58,17 @@ int nuc_cesam(const composition_map &comp,const matrix &T,const matrix &rho,
5858
int nuc_cesam_dcomp(composition_map &comp,const matrix &T,const matrix &rho,
5959
nuc_struct &nuc);
6060

61-
int eos_ideal(const matrix &X,double Z,const matrix &T,const matrix &p,
62-
matrix &rho,eos_struct &eos);
63-
int eos_idealrad(const matrix &X,double Z,const matrix &T,const matrix &p,
64-
matrix &rho,eos_struct &eos);
65-
int eos_opal(const matrix &X,double Z,const matrix &T,const matrix &p,
66-
matrix &rho,eos_struct &eos);
67-
int eos_freeeos(const matrix &X, double Z, const matrix &T, const matrix &p,
61+
int eos_ideal(const composition_map &, const matrix &T, const matrix &p,
62+
matrix &rho, eos_struct &eos);
63+
int eos_idealrad(const composition_map &, const matrix &T, const matrix &p,
64+
matrix &rho, eos_struct &eos);
65+
int eos_opal(const composition_map &, const matrix &T, const matrix &p,
66+
matrix &rho, eos_struct &eos);
67+
int eos_freeeos(const composition_map &, const matrix &T, const matrix &p,
6868
matrix &rho, eos_struct &eos);
6969

70-
int atm_onelayer(const matrix &X,double Z,const matrix &g,const matrix &Teff,
71-
const char *eos_name,const char *opa_name,atm_struct &atm);
70+
int atm_onelayer(const composition_map &, const matrix &g, const matrix &Teff,
71+
const char *eos_name, const char *opa_name, atm_struct &atm);
7272

7373
#endif
7474

‎src/physics/atm_onelayer.cpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -7,8 +7,8 @@
77
#include <string.h>
88
#include <cmath>
99

10-
int atm_onelayer(const matrix &X,double Z,const matrix &g,const matrix &Teff,
11-
const char *eos_name,const char *opa_name,atm_struct &atm) {
10+
int atm_onelayer(const composition_map &chemical_comp, const matrix &g, const matrix &Teff,
11+
const char *eos_name, const char *opa_name, atm_struct &atm) {
1212

1313
int n=g.nrows();
1414
int m=g.ncols();
@@ -43,8 +43,8 @@ int atm_onelayer(const matrix &X,double Z,const matrix &g,const matrix &Teff,
4343
}
4444
double F,dF,dlogps;
4545
p=pow(10,logps)*ones(1,1);
46-
if(eos_calc(X(i,j)*ones(1,1),Z,T,p,rho,eos)) return 1;
47-
if(opa_calc(X(i,j)*ones(1,1),Z,T,rho,opa)) return 1;
46+
if(eos_calc(chemical_comp(i,j)*ones(1,1), T, p, rho, eos)) return 1;
47+
if(opa_calc(chemical_comp.X()(i,j)*ones(1,1), chemical_comp.Z()(i,j), T, rho, opa)) return 1;
4848
F=logps+log10(opa.k)(0)-logg-log10(2./3.);
4949
dF=1.-(1+opa.dlnxi_lnrho(0))/eos.chi_rho(0);
5050
dlogps=-F/dF;

‎src/physics/eos_freeeos.cpp‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -47,7 +47,7 @@ void set_eps(double* eps, const double_map &local_chemical_mix) {
4747
eps[19] = local_chemical_mix["Ni"] / 58.693; // Ni //there are too many isotops
4848
}
4949

50-
int eos_freeeos(const matrix &X, double Z, const matrix &T, const matrix &p,
50+
int eos_freeeos(const composition_map &chemical_comp, const matrix &T, const matrix &p,
5151
matrix &rho, eos_struct &eos) {
5252

5353
double t;
@@ -127,12 +127,13 @@ int eos_freeeos(const matrix &X, double Z, const matrix &T, const matrix &p,
127127
&iteration_count);
128128
rho(i) = rhoi;
129129
if (iteration_count < 0) {
130+
// TODO: this makes no sense freeEOS has no tables...
130131
ester_err(
131132
"Values outside freeEOS eos table:\n"
132133
" X = %e\n"
133134
" Z = %e\n"
134135
" T = %e\n"
135-
" p = %e", X(i), Z, t, p(i));
136+
" p = %e", chemical_comp.X()(i), chemical_comp.Z()(i), t, p(i));
136137
}
137138
eos.s(i) = entropy[0];
138139
eos.G1(i) = gamma1;

‎src/physics/eos_ideal.cpp‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -5,12 +5,13 @@
55
#include "physics.h"
66
#include "constants.h"
77

8-
int eos_ideal(const matrix &X,double Z,const matrix &T,const matrix &p,
9-
matrix &rho,eos_struct &eos) {
8+
int eos_ideal(const composition_map &chemical_comp, const matrix &T, const matrix &p,
9+
matrix &rho, eos_struct &eos) {
1010

1111
matrix mu,b;
12-
13-
mu=4/(3+5*X-Z);
12+
13+
// TODO: this mix between X and Z implicitly assume a certain Z mixture
14+
mu = 4/(3 + 5*chemical_comp.X() - chemical_comp.Z());
1415
eos.prad=zeros(T.nrows(),T.ncols());
1516
rho=p/T/K_BOL*mu*HYDROGEN_MASS;
1617
b=ones(T.nrows(),T.ncols());

‎src/physics/eos_idealrad.cpp‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -4,12 +4,13 @@
44
#include "physics.h"
55
#include "constants.h"
66

7-
int eos_idealrad(const matrix &X,double Z,const matrix &T,const matrix &p,
8-
matrix &rho,eos_struct &eos) {
7+
int eos_idealrad(const composition_map &chemical_comp, const matrix &T, const matrix &p,
8+
matrix &rho, eos_struct &eos) {
99

1010
matrix mu,b;
1111

12-
mu=4/(3+5*X-Z);
12+
// TODO: this mix between X and Z implicitly assume a certain Z mixture
13+
mu = 4/(3 + 5*chemical_comp.X() - chemical_comp.Z());
1314
eos.prad=A_RAD/3*pow(T,4);
1415
rho=(p-eos.prad)/T/K_BOL*mu*HYDROGEN_MASS;
1516
b=(p-eos.prad)/p;

‎src/physics/eos_opal.cpp‎

Lines changed: 9 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -14,27 +14,16 @@ extern"C" {
1414
extern struct{
1515
double esact,eos[10];
1616
} eeos_;
17-
extern struct{
18-
int itime;
19-
} lreadco_;
2017
}
2118

22-
int eos_opal(const matrix &X,double Z,const matrix &T,const matrix &p,
23-
matrix &rho,eos_struct &eos) {
19+
int eos_opal(const composition_map &chemical_comp, const matrix &T, const matrix &p,
20+
matrix &rho, eos_struct &eos) {
2421

2522
matrix t6,p_mb;
2623
int i,N,error=0;
2724

28-
static double Z_table=-99;
29-
3025
t6=T*1e-6;
3126
p_mb=p*1e-12;
32-
33-
if(Z!=Z_table) {
34-
lreadco_.itime=0;
35-
zfs_interp_eos5_(&Z);
36-
Z_table=Z;
37-
}
3827

3928
if(error) {
4029
printf("Can't initialize OPAL EOS table\n");
@@ -58,8 +47,13 @@ int eos_opal(const matrix &X,double Z,const matrix &T,const matrix &p,
5847

5948
double Xi,Zi,t6i,p_mbi,rhoi;
6049
for(i=0;i<N;i++) {
61-
Xi=X(i);Zi=Z;t6i=t6(i);p_mbi=p_mb(i);
62-
eos5_xtrin_(&Xi,&Zi,&t6i,&p_mbi,&rhoi);
50+
Xi = chemical_comp.X()(i);
51+
Zi = chemical_comp.Z()(i);
52+
t6i = t6(i);
53+
p_mbi = p_mb(i);
54+
55+
zfs_interp_eos5_(&Zi);
56+
eos5_xtrin_(&Xi,&Zi,&t6i,&p_mbi,&rhoi);
6357
rho(i)=rhoi;
6458
eos.s(i)=1e6*(*(eeos_.eos+2));
6559
eos.G1(i)=*(eeos_.eos+7);

‎src/physics/physics.cpp‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -26,19 +26,19 @@ int opa_calc(const matrix &X,double Z,const matrix &T,const matrix &rho,
2626
return error;
2727
}
2828

29-
int eos_calc(const matrix &X,double Z,const matrix &T,const matrix &p,
30-
matrix &rho,eos_struct &eos) {
29+
int eos_calc(const composition_map &chemical_comp, const matrix &T, const matrix &p,
30+
matrix &rho, eos_struct &eos) {
3131

3232
int error=0;
3333

3434
if(!strcmp(eos.name,"ideal"))
35-
error=eos_ideal(X,Z,T,p,rho,eos);
35+
error = eos_ideal(chemical_comp, T, p, rho, eos);
3636
else if(!strcmp(eos.name,"ideal+rad"))
37-
error=eos_idealrad(X,Z,T,p,rho,eos);
37+
error = eos_idealrad(chemical_comp, T, p, rho, eos);
3838
else if(!strcmp(eos.name,"opal"))
39-
error=eos_opal(X,Z,T,p,rho,eos);
39+
error = eos_opal(chemical_comp, T, p, rho, eos);
4040
else if(!strcmp(eos.name,"freeeos"))
41-
error = eos_freeeos(X, Z, T, p, rho, eos);
41+
error = eos_freeeos(chemical_comp, T, p, rho, eos);
4242
else {
4343
ester_err("Unknown equation of state: %s",eos.name);
4444
return 1;
@@ -66,13 +66,13 @@ int nuc_calc(const composition_map &comp, const matrix &T, const matrix &rho,
6666

6767
}
6868

69-
int atm_calc(const matrix &X,double Z,const matrix &g,const matrix &Teff,
70-
const char *eos_name,const char *opa_name,atm_struct &atm) {
69+
int atm_calc(const composition_map &chemical_comp, const matrix &g, const matrix &Teff,
70+
const char *eos_name, const char *opa_name, atm_struct &atm) {
7171

7272
int error=0;
7373

7474
if(!strcmp(atm.name,"onelayer")) {
75-
error=atm_onelayer(X,Z,g,Teff,eos_name,opa_name,atm);
75+
error = atm_onelayer(chemical_comp, g, Teff, eos_name, opa_name, atm);
7676
} else {
7777
ester_err("Unknown atmosphere type: %s", atm.name);
7878
return 1;

‎src/star/star_phys.cpp‎

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -35,10 +35,12 @@ void star2d::eq_state() {
3535
int error;
3636

3737
matrix rhoc_m(1,1);
38-
eos_calc(comp.X()(0,0)*ones(1,1),Z0,ones(1,1)*Tc,ones(1,1)*pc,rhoc_m,eos);
39-
rhoc=rhoc_m(0);
38+
// TODO: is the following really usefull???
39+
// couldn't we just get rhoc from rho(0) ??? instead of this weird call only computing central value
40+
eos_calc(comp(0,0)*ones(1,1), Tc*ones(1,1), pc*ones(1,1), rhoc_m, eos);
41+
rhoc = rhoc_m(0);
4042

41-
error=eos_calc(comp.X(),Z0,T*Tc,p*pc,rho,eos);
43+
error = eos_calc(comp, T*Tc, p*pc, rho, eos);
4244

4345

4446
// rhoc=rho(0);
@@ -61,7 +63,7 @@ void star2d::atmosphere() {
6163
}
6264

6365

64-
error=atm_calc(comp.X(),Z0,gsup(),Teff(),eos.name,opa.name,atm);
66+
error = atm_calc(comp, gsup(), Teff(), eos.name, opa.name, atm);
6567

6668
if(error) exit(1);
6769

0 commit comments

Comments
 (0)