Skip to content

Commit

Permalink
background.c, chemequil.c, fpehandler.c, giigen.c, hydrogen.c, initia…
Browse files Browse the repository at this point in the history 
…l_xdr.c, iterate.c, kurucz.c,

ltepops.c, molzeeman.c, planck.c, readatom.c, readinput.c, readj.c, readmolecule.c, readvalue.c, 
scatter.c, atmos.h, atom.h, constant.h, inputs.h:
     Update to Han's latest RH version.
  • Loading branch information
@tiagopereira
tiagopereira committed Jan 31, 2012
1 parent 1444491 commit 0226ac0
Show file tree
Hide file tree
Showing 21 changed files with 820 additions and 776 deletions.
4 changes: 2 additions & 2 deletions atmos.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
Version: rh2.0
Author: Han Uitenbroek ([email protected])
Last modified: Wed Mar 3 10:56:48 2010 --
Last modified: Fri Jul 8 15:36:04 2011 --
-------------------------- ----------RH-- */

Expand Down Expand Up @@ -74,7 +74,7 @@ flags passive_bb(double lambda, int nspect, int mu, bool_t to_obs,
double *chi, double *eta, double *chip);

flags rlk_opacity(double lambda, int nspect, int mu, bool_t to_obs,
double *chi, double *eta, double *chip);
double *chi, double *eta, double *scatt, double *chip);

flags MolecularOpacity(double lambda, int nspect, int mu, bool_t to_obs,
double *chi, double *eta, double *chip);
Expand Down
5 changes: 3 additions & 2 deletions atom.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
Version: rh2.0
Author: Han Uitenbroek ([email protected])
Last modified: Tue Aug 4 16:26:49 2009 --
Last modified: Mon Apr 18 06:31:57 2011 --
-------------------------- ----------RH-- */

Expand Down Expand Up @@ -30,7 +30,7 @@ enum type {ATOMIC_LINE, ATOMIC_CONTINUUM,
enum ftype {FIXED_LINE, FIXED_CONTINUUM};
enum Trad_option {TRAD_ATMOSPHERIC, TRAD_PHOTOSPHERIC, TRAD_CHROMOSPHERIC};
enum vdWaals {UNSOLD, RIDDER_RENSBERGEN, BARKLEM, KURUCZ};
enum fit_type {KURUCZ_70, KURUCZ_85, SAUVAL_TATUM_84, TSUJI_73};
enum fit_type {KURUCZ_70, KURUCZ_85, SAUVAL_TATUM_84, IRWIN_81, TSUJI_73};
enum Hund {CASE_A, CASE_B};
enum Barklemtype {SP, PD, DF};
enum orbit_am {S_ORBIT=0, P_ORBIT, D_ORBIT, F_ORBIT};
Expand Down Expand Up @@ -119,6 +119,7 @@ struct Atom {
enum solution initial_solution;
int Nlevel, Nline, Ncont, Nfixed, Nprd, *stage, periodic_table,
activeindex;
long offset_coll;
double abundance, weight, *g, *E, **C, *vbroad, **n, **nstar,
*ntotal, **Gamma;
AtomicLine *line;
Expand Down
69 changes: 49 additions & 20 deletions background.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
Version: rh2.0
Author: Han Uitenbroek ([email protected])
Last modified: Fri Mar 5 09:16:05 2010 --
Last modified: Thu Sep 29 21:36:56 2011 --
-------------------------- ----------RH-- */

Expand Down Expand Up @@ -84,7 +84,7 @@
In static atmospheres these quantities are just mapped to
atmos.chi_c and atmos.eta_c to save memory space.
Note: If analyzeoutput == FALSE the auxiliary output files for
Note: If write_analyze_output == FALSE the auxiliary output files for
the Background Record Structure (BRS), metals, and molecules
are NOT written. This option is used when Background is called
from programs like solveray (formal solution along one specific
Expand Down Expand Up @@ -137,16 +137,16 @@ extern char messageStr[];

/* ------- begin -------------------------- Background.c ------------ */

void Background(bool_t analyzeoutput, bool_t equilibria_only)
void Background(bool_t write_analyze_output, bool_t equilibria_only)
{
const char routineName[] = "Background";
register int k, nspect, n, mu, to_obs;

static int ne_iter = 0;
char inputLine[MAX_LINE_SIZE];
bool_t exit_on_EOF, do_fudge, fromscratch;
int backgrrecno, index, Nfudge, NrecStokes;
double *chi, *chi_Q, *chi_U, *chi_V, *eta, *eta_Q, *eta_U, *eta_V,
*scatt, wavelength, *thomson, *chi_ai, *eta_ai,
double *chi, *eta, *scatt, wavelength, *thomson, *chi_ai, *eta_ai, *sca_ai,
Hmin_fudge, scatt_fudge, metal_fudge, *lambda_fudge, **fudge,
*Bnu, *chi_c, *eta_c, *sca_c, *chip, *chip_c;
Atom *He;
Expand All @@ -155,7 +155,13 @@ void Background(bool_t analyzeoutput, bool_t equilibria_only)

getCPU(2, TIME_START, NULL);

if (input.solve_ne) Solve_ne(atmos.ne, fromscratch=TRUE);
if (input.solve_ne == ONCE || input.solve_ne == ITERATION ) {
fromscratch = (input.solve_ne == ONCE ||
(input.solve_ne == ITERATION && ne_iter == 0)) ?
TRUE : FALSE;
Solve_ne(atmos.ne, fromscratch);
ne_iter++;
}
SetLTEQuantities();

if (input.NonICE)
Expand Down Expand Up @@ -249,9 +255,11 @@ void Background(bool_t analyzeoutput, bool_t equilibria_only)
if (atmos.moving || atmos.Stokes) {
chi_ai = (double *) malloc(atmos.Nspace * sizeof(double));
eta_ai = (double *) malloc(atmos.Nspace * sizeof(double));
sca_ai = (double *) malloc(atmos.Nspace * sizeof(double));
} else {
chi_ai = chi_c;
eta_ai = eta_c;
sca_ai = sca_c;
}
Bnu = (double *) malloc(atmos.Nspace * sizeof(double));

Expand Down Expand Up @@ -289,9 +297,9 @@ void Background(bool_t analyzeoutput, bool_t equilibria_only)

if (atmos.moving || atmos.Stokes) {
atmos.backgrrecno =
(long *) malloc(2*spectrum.Nspect*atmos.Nrays * sizeof(long));
(long *) malloc(2*spectrum.Nspect*atmos.Nrays * sizeof(int));
} else
atmos.backgrrecno = (long *) malloc(spectrum.Nspect * sizeof(long));
atmos.backgrrecno = (long *) malloc(spectrum.Nspect * sizeof(int));

/* --- Open output file for background opacity, emissivity,
scattering -- -------------- */
Expand Down Expand Up @@ -323,7 +331,7 @@ void Background(bool_t analyzeoutput, bool_t equilibria_only)
for (k = 0; k < atmos.Nspace; k++) {
chi_ai[k] = 0.0;
eta_ai[k] = 0.0;
sca_c[k] = thomson[k];
sca_ai[k] = thomson[k];
}
/* --- Negative hydrogen ion, bound-free and free-free -- ------- */

Expand Down Expand Up @@ -380,14 +388,14 @@ void Background(bool_t analyzeoutput, bool_t equilibria_only)

if (Rayleigh(wavelength, atmos.H, scatt)) {
for (k = 0; k < atmos.Nspace; k++) {
sca_c[k] += scatt[k];
sca_ai[k] += scatt[k];
}
}
/* --- Rayleigh scattering by neutral helium -- -------------- */

if (He && Rayleigh(wavelength, He, scatt)) {
for (k = 0; k < atmos.Nspace; k++) {
sca_c[k] += scatt[k];
sca_ai[k] += scatt[k];
}
}
/* --- Absorption by H + H^+ (referred to as H2plus free-free) -- */
Expand All @@ -403,7 +411,7 @@ void Background(bool_t analyzeoutput, bool_t equilibria_only)

if (Rayleigh_H2(wavelength, scatt)) {
for (k = 0; k < atmos.Nspace; k++) {
sca_c[k] += scatt[k];
sca_ai[k] += scatt[k];
}
}
if (H2minus_ff(wavelength, chi)) {
Expand Down Expand Up @@ -437,9 +445,10 @@ void Background(bool_t analyzeoutput, bool_t equilibria_only)
} else {
scatt_fudge = 1.0;
}
for (k = 0; k < atmos.Nspace; k++)
chi_ai[k] += sca_c[k] * scatt_fudge;

for (k = 0; k < atmos.Nspace; k++) {
sca_ai[k] *= scatt_fudge;
chi_ai[k] += sca_ai[k];
}
/* --- Now the contributions that may be angle-dependent due to the
presence of atomic or molecular lines -- -------------- */

Expand All @@ -453,7 +462,9 @@ void Background(bool_t analyzeoutput, bool_t equilibria_only)
for (k = 0; k < atmos.Nspace; k++) {
chi_c[k] = chi_ai[k];
eta_c[k] = eta_ai[k];
sca_c[k] = sca_ai[k];
}

/* --- Zero the polarized quantities, if necessary -- ----- */

if (atmos.Stokes) {
Expand Down Expand Up @@ -493,7 +504,7 @@ void Background(bool_t analyzeoutput, bool_t equilibria_only)

if (atmos.Nrlk > 0) {
backgrflags = rlk_opacity(wavelength, nspect, mu, to_obs,
chi, eta, chip);
chi, eta, scatt, chip);
if (backgrflags.hasline) {
atmos.backgrflags[nspect].hasline = TRUE;
if (backgrflags.ispolarized) {
Expand All @@ -510,6 +521,12 @@ void Background(bool_t analyzeoutput, bool_t equilibria_only)
chi_c[k] += chi[k];
eta_c[k] += eta[k];
}
if (input.rlkscatter) {
for (k = 0; k < atmos.Nspace; k++) {
sca_c[k] += scatt[k];
chi_c[k] += scatt[k];
}
}
}
}
/* --- Add opacity from molecular lines -- -------------- */
Expand Down Expand Up @@ -564,13 +581,19 @@ void Background(bool_t analyzeoutput, bool_t equilibria_only)

if (atmos.Nrlk > 0) {
backgrflags = rlk_opacity(wavelength, nspect, 0, TRUE,
chi, eta, NULL);
chi, eta, scatt, NULL);
if (backgrflags.hasline) {
atmos.backgrflags[nspect].hasline = TRUE;
for (k = 0; k < atmos.Nspace; k++) {
chi_c[k] += chi[k];
eta_c[k] += eta[k];
}
if (input.rlkscatter) {
for (k = 0; k < atmos.Nspace; k++) {
sca_c[k] += scatt[k];
chi_c[k] += scatt[k];
}
}
}
}
/* --- Add opacity from molecular lines -- -------------- */
Expand All @@ -592,7 +615,7 @@ void Background(bool_t analyzeoutput, bool_t equilibria_only)
}
}

if (analyzeoutput) {
if (write_analyze_output) {
/* --- Write background record structure -- ------------ */

writeBRS();
Expand All @@ -607,12 +630,17 @@ void Background(bool_t analyzeoutput, bool_t equilibria_only)

if (atmos.Natom > 1) {
for (n = 1; n < atmos.Natom; n++)
if (!atmos.atoms[n].active && !atmos.hydrostatic)
if (!atmos.atoms[n].active &&
!atmos.hydrostatic &&
input.solve_ne < ITERATION)
freeAtom(&atmos.atoms[n]);
}
if (atmos.Nmolecule > 1) {
for (n = 1; n < atmos.Nmolecule; n++)
if (!atmos.molecules[n].active) freeMolecule(&atmos.molecules[n]);
if (!atmos.molecules[n].active &&
!atmos.hydrostatic &&
input.solve_ne < ITERATION)
freeMolecule(&atmos.molecules[n]);
}

if (strcmp(input.KuruczData, "none")) {
Expand All @@ -638,6 +666,7 @@ void Background(bool_t analyzeoutput, bool_t equilibria_only)
if (atmos.moving || atmos.Stokes) {
free(chi_ai);
free(eta_ai);
free(sca_ai);
}
if (atmos.Stokes && input.magneto_optical) {
free(chip);
Expand Down
20 changes: 15 additions & 5 deletions chemequil.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
Version: rh2.0
Author: Han Uitenbroek ([email protected])
Last modified: Tue Mar 31 09:15:25 2009 --
Last modified: Mon Apr 18 07:08:41 2011 --
-------------------------- ----------RH-- */

Expand All @@ -23,9 +23,6 @@

#define COMMENT_CHAR "#"

/* --- Ionization energy Hmin in [J] -- -------------- */

#define E_ION_HMIN 0.754 * EV

/* --- Acceleration parameters -- -------------- */

Expand Down Expand Up @@ -342,7 +339,7 @@ void ChemicalEquilibrium(int NmaxIter, double iterLimit)
}
/* --- Store Hmin density -- -------------- */

atmos.nHmin[k] = atmos.ne[k] * n[0] * PhiHmin;
atmos.nHmin[k] = atmos.ne[k] * (n[0] * PhiHmin);

/* --- Store molecular densities -- -------------- */

Expand Down Expand Up @@ -428,6 +425,15 @@ double partfunction(struct Molecule *molecule, double T)
pf = POW10(pf);
break;

case IRWIN_81:
t = log(T);
pf = molecule->pf_coef[0];
for (i = 1; i < molecule->Npf; i++)
pf = pf*t + molecule->pf_coef[i];

pf = exp(pf);
break;

case TSUJI_73:
break;
default:
Expand Down Expand Up @@ -481,6 +487,10 @@ double equilconstant(struct Molecule *molecule, double T)
cgs_to_SI = pow(CUBE(CM_TO_M), mk);
break;

case IRWIN_81:

/* ---- Use SAUVAL_TATUM_84 for equilibrium constant -- ---------- */
;
case SAUVAL_TATUM_84:
theta = THETA0 / T;
t = log10(theta);
Expand Down
7 changes: 6 additions & 1 deletion constant.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
Version: rh2.0
Author: Han Uitenbroek ([email protected])
Last modified: Wed Sep 14 22:13:01 2005 --
Last modified: Wed Nov 17 16:27:27 2010 --
-------------------------- ----------RH-- */

Expand Down Expand Up @@ -60,13 +60,18 @@
#endif
#define SQRTPI 1.77245385090551


/* --- 1/(2sqrt(2)), needed for anisotropy of radiation -- ---------- */

#define TWOSQRTTWO 0.35355339059327

#define LARMOR (Q_ELECTRON / (4.0*PI*M_ELECTRON)) * NM_TO_M


/* --- Ionization energy Hmin in [J] -- -------------- */

#define E_ION_HMIN 0.754 * EV

#endif /* !__CONSTANT_H__ */

/* ------- end ---------------------------- constant.h -------------- */
Loading

0 comments on commit 0226ac0

Please sign in to comment.