Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix Barklem broadening for Kurucz lines, add support for JJ and JK coupling #35

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Fix Barklem broadening for Kurucz lines, add support for JJ and JK coupling #35

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8dcac41
Added updates from Han for Zeeman calculations
tiagopereira Nov 11, 2021
b041a88
Fix for correct Barklem broadening for Kurucz lines, general cleanup
tiagopereira Nov 11, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 22 additions & 11 deletions atom.h
View file
Open in desktop
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: Mon Apr 18 06:31:57 2011 --
Last modified: Tue May 25 18:13:07 2021 --

-------------------------- ----------RH-- */

Expand All @@ -22,8 +22,6 @@
#define PRD_QSPREAD 5.0
#define PRD_DQ 0.25

#define RLK_LABEL_LENGTH 10


enum type {ATOMIC_LINE, ATOMIC_CONTINUUM,
VIBRATION_ROTATION, MOLECULAR_ELECTRONIC};
Expand All @@ -34,6 +32,7 @@ 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};
enum zeeman_cpl {LS_COUPLING=0, JK_COUPLING, JJ_COUPLING};

/* --- Structure prototypes -- -------------- */

Expand Down Expand Up @@ -63,6 +62,7 @@ struct AtomicLine {
Atom *atom;
AtomicLine **xrd;
pthread_mutex_t rate_lock;
ZeemanMultiplet *zm;
};

typedef struct {
Expand Down Expand Up @@ -157,20 +157,27 @@ struct Molecule {
pthread_mutex_t Gamma_lock;
};

typedef struct {
int L, L1, l1, l2, l, L_ac;
double g, E, S, J, S1, J1, j1, j2, K, gL, hfs;
enum zeeman_cpl cpl;
bool_t zm_explicit;
} RLK_level;

typedef struct {
bool_t polarizable;
enum vdWaals vdwaals;
int pt_index, stage, isotope, Li, Lj;
double lambda0, gi, gj, Ei, Ej, Bji, Aji, Bij, Si, Sj,
int pt_index, stage, isotope;
double lambda0, Bji, Aji, Bij,
Grad, GStark, GvdWaals, hyperfine_frac,
isotope_frac, gL_i, gL_j, hfs_i, hfs_j, iso_dl,
cross, alpha;
isotope_frac, iso_dl, cross, alpha;
RLK_level level_i, level_j;
ZeemanMultiplet *zm;
} RLK_Line;

struct ZeemanMultiplet{
int Ncomponent, *q;
double *shift, *strength;
double *shift, *strength, g_eff;
};

typedef struct {
Expand Down Expand Up @@ -211,6 +218,10 @@ void initAtom(Atom *atom);
void initAtomicLine(AtomicLine *line);
void initAtomicContinuum(AtomicContinuum *continuum);

void initZeeman(ZeemanMultiplet *zm);
void freeZeeman(ZeemanMultiplet *zm);
double zm_gamma(double J, double S, double L);

void initGammaAtom(Atom *atom, double cswitch);
void initGammaMolecule(Molecule *molecule);

Expand Down Expand Up @@ -287,9 +298,9 @@ bool_t determinate(char *label, double g, int *n, double *S, int *L,
double effectiveLande(AtomicLine *line);
double Lande(double S, int L, double J);
void StokesProfile(AtomicLine *line);
ZeemanMultiplet* Zeeman(AtomicLine *line);
ZeemanMultiplet* MolZeeman(MolecularLine *mrt);
double MolLande_eff(MolecularLine *mrt);
void Zeeman(AtomicLine *line);
void MolZeeman(MolecularLine *mrt);
double MolLande_eff(MolecularLine *mrt);
int getOrbital(char orbit);
double ZeemanStrength(double Ju, double Mu, double Jl, double Ml);

Expand Down
6 changes: 4 additions & 2 deletions background.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -190,6 +190,9 @@ void Background(bool_t write_analyze_output, bool_t equilibria_only)
}

getCPU(3, TIME_START, NULL);

do_fudge = FALSE;

if (strcmp(input.fudgeData, "none")) {
do_fudge = TRUE;

Expand All @@ -216,8 +219,7 @@ void Background(bool_t write_analyze_output, bool_t equilibria_only)
}
for (n = 0; n < 3*Nfudge; n++) fudge[0][n] += 1.0;
fclose(fp_fudge);
} else
do_fudge = FALSE;
}

/* --- Allocate temporary storage space. The quantities are used
for the following purposes:
Expand Down
11 changes: 5 additions & 6 deletions barklem.c
View file
Open in desktop
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 Nov 5, 2014, 15:24 --
Last modified: Fri May 14 19:18:30 2021 --

-------------------------- ----------RH-- */

Expand Down Expand Up @@ -173,18 +173,18 @@ bool_t getBarklemcross(Barklemstruct *bs, RLK_Line *rlk)
if (rlk->stage > 0)
return FALSE;

if ((deltaEi = element->ionpot[rlk->stage] - rlk->Ei) <= 0.0)
if ((deltaEi = element->ionpot[rlk->stage] - rlk->level_i.E) <= 0.0)
return FALSE;
if ((deltaEj = element->ionpot[rlk->stage] - rlk->Ej) <= 0.0)
if ((deltaEj = element->ionpot[rlk->stage] - rlk->level_j.E) <= 0.0)
return FALSE;

Z = (double) (rlk->stage + 1);
E_Rydberg = E_RYDBERG / (1.0 + M_ELECTRON / (element->weight * AMU));
neff1 = Z * sqrt(E_Rydberg / deltaEi);
neff2 = Z * sqrt(E_Rydberg / deltaEj);

if (rlk->Li > rlk->Lj) SWAPDOUBLE(neff1, neff2);

if (rlk->level_i.L_ac > rlk->level_j.L_ac) SWAPDOUBLE(neff1, neff2);
if (neff1 < bs->neff1[0] || neff1 > bs->neff1[bs->N1-1])
return FALSE;
Locate(bs->N1, bs->neff1, neff1, &index);
Expand All @@ -204,7 +204,6 @@ bool_t getBarklemcross(Barklemstruct *bs, RLK_Line *rlk)
rlk->alpha = cubeconvol(bs->N2, bs->N1,
bs->alpha[0], findex2, findex1);


reducedmass = AMU / (1.0/atmos.H->weight + 1.0/element->weight);
meanvelocity = sqrt(8.0 * KBOLTZMANN / (PI * reducedmass));
crossmean = SQ(RBOHR) * pow(meanvelocity / 1.0E4, -rlk->alpha);
Expand Down
2 changes: 0 additions & 2 deletions collision.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -957,8 +957,6 @@ void CollisionRate(struct Atom *atom, char *fp_atom)
free(T);
free(coeff);

//fsetpos(fp_atom, &collpos); // Tiago: not needed now (but when to free fp_atom?!)

sprintf(labelStr, "Collision Rate %2s", atom->ID);
getCPU(3, TIME_POLL, labelStr);
}
Expand Down
4 changes: 2 additions & 2 deletions fpehandler.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -105,8 +105,8 @@ void SetFPEtraps(void)
/* --- Enable some exceptions.
At startup all exceptions are masked. -- -------------- */

/* Tiago: commented this out for the 1.5D version (where it should
not stop execution!)
/* Commented out for the 1.5D version, so it does not stop execution.

feenableexcept(FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW); */
}

Expand Down
2 changes: 1 addition & 1 deletion inputs.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@ typedef struct {
bool_t magneto_optical, XRD, Eddington,
backgr_pol, limit_memory, allow_passive_bb, NonICE,
rlkscatter, prdh_limit_mem, backgr_in_mem, xdr_endian,
old_background, accelerate_mols;
old_background, accelerate_mols, RLK_explicit;
enum solution startJ;
enum StokesMode StokesMode;
enum S_interpol S_interpolation;
Expand Down
42 changes: 4 additions & 38 deletions iterate.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,8 +56,6 @@ void Iterate(int NmaxIter, double iterLimit)
double dpopsmax, PRDiterlimit;
Atom *atom;
Molecule *molecule;
AtomicLine *line; // Tiago: DELETE
int i, mu, to_obs, lamu; // Tiago: DELETE

if (NmaxIter <= 0) return;
getCPU(1, TIME_START, NULL);
Expand Down Expand Up @@ -119,9 +117,9 @@ void Iterate(int NmaxIter, double iterLimit)
/* --- Redistribute intensity in PRD lines if necessary -- ---- */

if (input.PRDiterLimit < 0.0)
PRDiterlimit = MAX(dpopsmax, -input.PRDiterLimit);
PRDiterlimit = MAX(dpopsmax, -input.PRDiterLimit);
else
PRDiterlimit = input.PRDiterLimit;
PRDiterlimit = input.PRDiterLimit;
Redistribute(input.PRD_NmaxIter, PRDiterlimit);
}

Expand All @@ -146,46 +144,14 @@ void Iterate(int NmaxIter, double iterLimit)

if (atmos.hydrostatic) {
if (!atmos.atoms[0].active) {
sprintf(messageStr, "Can only perform hydrostatic equilibrium"
sprintf(messageStr, "Can only perform hydrostatic equilibrium"
" for hydrogen active");
Error(ERROR_LEVEL_2, routineName, messageStr);
Error(ERROR_LEVEL_2, routineName, messageStr);
}
Hydrostatic(N_MAX_HSE_ITER, HSE_ITER_LIMIT);
}
}

// Tiago: temporary printouts to get PRD rho after iteration
/*
atom = atmos.activeatoms[0];
line = &atom->line[0];

switch (input.PRD_angle_dep) {
case PRD_ANGLE_INDEP:
printf("rho_prd = \n");
for (i = 0; i < line->Nlambda; i++) {
printf("%8.4f %e %e %e %e %e\n", line->lambda[i], line->rho_prd[i][105], line->rho_prd[i][110], line->rho_prd[i][120], line->rho_prd[i][150], line->rho_prd[i][155]);
}
//exit(1);
break;

case PRD_ANGLE_DEP:
for (mu = 0; mu < atmos.Nrays; mu++) {
for (to_obs = 0; to_obs <= 1; to_obs++) {
for (i = 0; i < line->Nlambda; i++) {
lamu = 2*(atmos.Nrays*i + mu) + to_obs;
if ((to_obs == 1) && (mu == 4))
printf("%8.4f %e %e %e %e %e\n", line->lambda[i], line->rho_prd[lamu][105],line->rho_prd[lamu][110],line->rho_prd[lamu][120], line->rho_prd[lamu][150], line->rho_prd[lamu][155] );
}
}
}
//exit(1);
break;
}
*/




for (nact = 0; nact < atmos.Nactiveatom; nact++) {
atom = atmos.activeatoms[nact];
freeMatrix((void **) atom->Gamma);
Expand Down
Loading