diff --git a/atom.h b/atom.h index 17e7458..2006751 100644 --- a/atom.h +++ b/atom.h @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Mon Apr 18 06:31:57 2011 -- + Last modified: Tue May 25 18:13:07 2021 -- -------------------------- ----------RH-- */ @@ -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}; @@ -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 -- -------------- */ @@ -63,6 +62,7 @@ struct AtomicLine { Atom *atom; AtomicLine **xrd; pthread_mutex_t rate_lock; + ZeemanMultiplet *zm; }; typedef struct { @@ -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 { @@ -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); @@ -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); diff --git a/background.c b/background.c index fe89f7b..3a56ab0 100644 --- a/background.c +++ b/background.c @@ -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; @@ -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: diff --git a/barklem.c b/barklem.c index 1680aa8..328c555 100644 --- a/barklem.c +++ b/barklem.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Wed Nov 5, 2014, 15:24 -- + Last modified: Fri May 14 19:18:30 2021 -- -------------------------- ----------RH-- */ @@ -173,9 +173,9 @@ 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); @@ -183,8 +183,8 @@ bool_t getBarklemcross(Barklemstruct *bs, RLK_Line *rlk) 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); @@ -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); diff --git a/collision.c b/collision.c index c94637d..e556d14 100644 --- a/collision.c +++ b/collision.c @@ -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); } diff --git a/fpehandler.c b/fpehandler.c index 00dc472..f3a9465 100644 --- a/fpehandler.c +++ b/fpehandler.c @@ -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); */ } diff --git a/inputs.h b/inputs.h index 111fe1f..cffeab9 100644 --- a/inputs.h +++ b/inputs.h @@ -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; diff --git a/iterate.c b/iterate.c index 5c2c311..4605a2d 100644 --- a/iterate.c +++ b/iterate.c @@ -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); @@ -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); } @@ -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); diff --git a/kurucz.c b/kurucz.c index d79bf8c..66cb4e6 100644 --- a/kurucz.c +++ b/kurucz.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Thu Oct 12 14:52:27 2017 -- + Last modified: Fri May 28 16:02:49 2021 -- -------------------------- ----------RH-- */ @@ -82,18 +82,23 @@ FORMAT(F11.4,F7.3,F6.2,F12.3,F5.2,1X,A10,F12.3,F5.2,1X,A10, #define MAX_GAUSS_DOPPLER 7.0 #define USE_TABULATED_WAVELENGTH 1 +#define RLK_LABEL_LENGTH 10 + /* --- Function prototypes -- -------------- */ -double RLKProfile(RLK_Line *rlk, int k, int mu, bool_t to_obs, - double lambda, - double *phi_Q, double *phi_U, double *phi_V, - double *psi_Q, double *psi_U, double *psi_V); -ZeemanMultiplet* RLKZeeman(RLK_Line *rlk); -void initRLK(RLK_Line *rlk); -bool_t RLKdeterminate(char *labeli, char *labelj, RLK_Line *rlk); -void getUnsoldcross(RLK_Line *rlk); -void free_BS(Barklemstruct *bs); +double RLKProfile(RLK_Line *rlk, int k, int mu, bool_t to_obs, + double lambda, + double *phi_Q, double *phi_U, double *phi_V, + double *psi_Q, double *psi_U, double *psi_V); +void RLKZeeman(RLK_Line *rlk); +double RLKLande(RLK_level* level); +void initRLK(RLK_Line *rlk); +bool_t RLKdet_level(char* label, RLK_level *level); +bool_t RLKdet_level_ac(char* label, RLK_level *level); +double getJK_K(char c); +void getUnsoldcross(RLK_Line *rlk); +void free_BS(Barklemstruct *bs); /* --- Global variables -- -------------- */ @@ -112,14 +117,14 @@ void readKuruczLines(char *inputFile) (Q_ELECTRON/M_ELECTRON) / CLIGHT; char inputLine[RLK_RECORD_LENGTH+1], listName[MAX_LINE_SIZE], - filename[MAX_LINE_SIZE], Gvalues[18+1], elem_code[7], + filename[MAX_LINE_SIZE], Gvalues[18+1], elem_code[7], labeli[RLK_LABEL_LENGTH+1], labelj[RLK_LABEL_LENGTH+1], *commentChar = COMMENT_CHAR; bool_t swap_levels, determined, useBarklem, exit_on_EOF; int Nline, Nread, Nrequired, checkPoint, hfs_i, hfs_j, gL_i, gL_j, - iso_dl, i; + iso_dl, Li, Lj, i; unsigned int line_index; - double lambda0, Ji, Jj, Grad, GStark, GvdWaals, pti, + double lambda0, Ji, Jj, Grad, GStark, GvdWaals, Ei, Ej, gf, lambda_air; size_t length; RLK_Line *rlk; @@ -139,18 +144,18 @@ void readKuruczLines(char *inputFile) /* Use saved or read file */ if (input.kurucz_file_contents != NULL) { - fp_Kurucz = input.kurucz_file_contents; + fp_Kurucz = input.kurucz_file_contents; } else { - fp_Kurucz = readWholeFile(inputFile); + fp_Kurucz = readWholeFile(inputFile); } /* Count line files, save content */ Nline = 0; file_string = fp_Kurucz; /* To avoid rewinding fp_Kurucz */ while (getLineString(&file_string, commentChar, listName, FALSE) != EOF) { - length = strlen(listName); - /* To avoid untraced bug that reads extra lines where there's only one */ - if (listName[length - 1] != '\n') continue; - Nline++; + length = strlen(listName); + /* To avoid untraced bug that reads extra lines where there's only one */ + if (listName[length - 1] != '\n') continue; + Nline++; } if (input.kurucz_file_contents == NULL) { /* Make space for input data */ input.kurucz_file_contents = fp_Kurucz; @@ -186,7 +191,7 @@ void readKuruczLines(char *inputFile) if (atmos.Nrlk == 0) atmos.rlk_lines = NULL; atmos.rlk_lines = (RLK_Line *) - realloc(atmos.rlk_lines, (Nline + atmos.Nrlk) * sizeof(RLK_Line)); + realloc(atmos.rlk_lines, (Nline + atmos.Nrlk) * sizeof(RLK_Line)); /* --- Read lines from file -- -------------- */ @@ -196,168 +201,155 @@ void readKuruczLines(char *inputFile) initRLK(rlk); - Nread = sscanf(inputLine, "%lf %lf %s %lf", - &lambda_air, &gf, (char *) &elem_code, &Ei); + Nread = sscanf(inputLine, "%lf %lf %s %lf", + &lambda_air, &gf, (char *) &elem_code, &Ei); /* --- Ionization stage and periodic table index -- --------- */ sscanf(elem_code, "%d.%d", &rlk->pt_index, &rlk->stage); - Nread += sscanf(inputLine+53, "%lf", &Ej); + Nread += sscanf(inputLine+53, "%lf", &Ej); - Ei = fabs(Ei) * (HPLANCK * CLIGHT) / CM_TO_M; - Ej = fabs(Ej) * (HPLANCK * CLIGHT) / CM_TO_M; + Ei = fabs(Ei) * (HPLANCK * CLIGHT) / CM_TO_M; + Ej = fabs(Ej) * (HPLANCK * CLIGHT) / CM_TO_M; - /* --- Beware: the Kurucz linelist has upper and lower levels - of a transition in random order. Therefore, we have to + /* --- Beware: the Kurucz linelist has upper and lower levels + of a transition in random order. Therefore, we have to check for the lowest energy of the two and use that as lower level -- -------------- */ - if (Ej < Ei) { - swap_levels = TRUE; - rlk->Ei = Ej; - rlk->Ej = Ei; - strncpy(labeli, inputLine+69, RLK_LABEL_LENGTH); - strncpy(labelj, inputLine+41, RLK_LABEL_LENGTH); - } else { - swap_levels = FALSE; - rlk->Ei = Ei; - rlk->Ej = Ej; - strncpy(labeli, inputLine+41, RLK_LABEL_LENGTH); - strncpy(labelj, inputLine+69, RLK_LABEL_LENGTH); - } - - Nread += sscanf(inputLine+35, "%lf", &Ji); - Nread += sscanf(inputLine+63, "%lf", &Jj); - if (swap_levels) SWAPDOUBLE(Ji, Jj); - rlk->gi = 2*Ji + 1; - rlk->gj = 2*Jj + 1; - - if (USE_TABULATED_WAVELENGTH) { - /* --- In this case use tabulated wavelength and adjust - upper-level energy -- -------------- */ - - air_to_vacuum(1, &lambda_air, &lambda0); - lambda0 *= NM_TO_M; - rlk->Ej = rlk->Ei + (HPLANCK * CLIGHT) / lambda0; - } else { - /* --- Else use energy levels to calculate lambda0 -- ----- */ - - lambda0 = (HPLANCK * CLIGHT) / (rlk->Ej - rlk->Ei); - } - rlk->Aji = C / SQ(lambda0) * POW10(gf) / rlk->gj; - rlk->Bji = CUBE(lambda0) / (2.0 * HPLANCK * CLIGHT) * rlk->Aji; - rlk->Bij = (rlk->gj / rlk->gi) * rlk->Bji; + if (Ej < Ei) { + swap_levels = TRUE; + rlk->level_i.E = Ej; + rlk->level_j.E = Ei; + strncpy(labeli, inputLine+69, RLK_LABEL_LENGTH); + strncpy(labelj, inputLine+41, RLK_LABEL_LENGTH); + } else { + swap_levels = FALSE; + rlk->level_i.E = Ei; + rlk->level_j.E = Ej; + strncpy(labeli, inputLine+41, RLK_LABEL_LENGTH); + strncpy(labelj, inputLine+69, RLK_LABEL_LENGTH); + } + + Nread += sscanf(inputLine+35, "%lf", &Ji); + Nread += sscanf(inputLine+63, "%lf", &Jj); + if (swap_levels) SWAPDOUBLE(Ji, Jj); + + rlk->level_i.J = Ji; + rlk->level_i.g = 2*Ji + 1; + rlk->level_j.J = Jj; + rlk->level_j.g = 2*Jj + 1; + + if (USE_TABULATED_WAVELENGTH) { + + /* --- In this case use tabulated wavelength and adjust + upper-level energy -- -------------- */ + + air_to_vacuum(1, &lambda_air, &lambda0); + lambda0 *= NM_TO_M; + rlk->level_j.E = rlk->level_i.E + (HPLANCK * CLIGHT) / lambda0; + } else { + /* --- Else use energy levels to calculate lambda0 -- ----- */ + + lambda0 = (HPLANCK * CLIGHT) / (rlk->level_j.E - rlk->level_i.E); + } + rlk->Aji = C / SQ(lambda0) * POW10(gf) / rlk->level_j.g; + rlk->Bji = CUBE(lambda0) / (2.0 * HPLANCK * CLIGHT) * rlk->Aji; + rlk->Bij = (rlk->level_j.g / rlk->level_i.g) * rlk->Bji; /* --- Store in nm -- -------------- */ + rlk->lambda0 = lambda0 / NM_TO_M; - rlk->lambda0 = lambda0 / NM_TO_M; - - /* --- Get quantum numbers for angular momentum and spin -- - */ - - determined = RLKdeterminate(labeli, labelj, rlk); + determined = (RLKdet_level(labeli, &rlk->level_i) && + RLKdet_level(labelj, &rlk->level_j)); rlk->polarizable = (atmos.Stokes && determined); + /* --- Get orbital quantum level for Barklem tables --- */ + determined = (RLKdet_level_ac(labeli, &rlk->level_i) && + RLKdet_level_ac(labelj, &rlk->level_j)); + /* --- Line broadening -- -------------- */ - strncpy(Gvalues, inputLine+79, 18); - Nread += sscanf(Gvalues, "%lf %lf %lf", &Grad, &GStark, &GvdWaals); + strncpy(Gvalues, inputLine+79, 18); + Nread += sscanf(Gvalues, "%lf %lf %lf", &Grad, &GStark, &GvdWaals); - if (GStark != 0.0) - rlk->GStark = POW10(GStark) * CUBE(CM_TO_M); - else - rlk->GStark = 0.0; + if (GStark != 0.0) + rlk->GStark = POW10(GStark) * CUBE(CM_TO_M); + else + rlk->GStark = 0.0; - if (GvdWaals != 0.0) - rlk->GvdWaals = POW10(GvdWaals) * CUBE(CM_TO_M); - else - rlk->GvdWaals = 0.0; + if (GvdWaals != 0.0) + rlk->GvdWaals = POW10(GvdWaals) * CUBE(CM_TO_M); + else + rlk->GvdWaals = 0.0; /* --- If possible use Barklem formalism -- -------------- */ - useBarklem = FALSE; - if (determined) { - if ((rlk->Li == S_ORBIT && rlk->Lj == P_ORBIT) || - (rlk->Li == P_ORBIT && rlk->Lj == S_ORBIT)) { - useBarklem = getBarklemcross(&bs_SP, rlk); - } else if ((rlk->Li == P_ORBIT && rlk->Lj == D_ORBIT) || - (rlk->Li == D_ORBIT && rlk->Lj == P_ORBIT)) { - useBarklem = getBarklemcross(&bs_PD, rlk); - } else if ((rlk->Li == D_ORBIT && rlk->Lj == F_ORBIT) || - (rlk->Li == F_ORBIT && rlk->Lj == D_ORBIT)) { - useBarklem = getBarklemcross(&bs_DF, rlk); - } - } - /* --- Else use good old Unsoeld -- -------------- */ + useBarklem = FALSE; + if (determined && + rlk->level_i.cpl == LS_COUPLING && + rlk->level_j.cpl == LS_COUPLING) { + + Li = rlk->level_i.L_ac; + Lj = rlk->level_j.L_ac; + + if ((Li == S_ORBIT && Lj == P_ORBIT) || + (Li == P_ORBIT && Lj == S_ORBIT)) { + useBarklem = getBarklemcross(&bs_SP, rlk); + } else if ((Li == P_ORBIT && Lj == D_ORBIT) || + (Li == D_ORBIT && Lj == P_ORBIT)) { + useBarklem = getBarklemcross(&bs_PD, rlk); + } else if ((Li == D_ORBIT && Lj == F_ORBIT) || + (Li == F_ORBIT && Lj == D_ORBIT)) { + useBarklem = getBarklemcross(&bs_DF, rlk); + } + } + /* --- Else use good old Unsoeld -- -------------- */ if (!useBarklem) { - getUnsoldcross(rlk); - } - /* --- Radiative broadening -- -------------- */ - - if (Grad != 0.0) { - rlk->Grad = POW10(Grad); - } else { - - /* --- Just take the Einstein Aji value, but only if either - Stark or vd Waals broadening is in effect -- ------- */ - - if (GStark != 0.0 || GvdWaals != 0.0) - rlk->Grad = rlk->Aji; - else { - rlk->Grad = 0.0; - - /* --- In this case the line is not polarizable because - there is no way to determine its damping -- ------ */ - - rlk->polarizable = FALSE; - } - } - /* --- Isotope and hyperfine fractions and slpittings -- ---- */ - - Nread += sscanf(inputLine+106, "%d", &rlk->isotope); - Nread += sscanf(inputLine+108, "%lf", &rlk->isotope_frac); - rlk->isotope_frac = POW10(rlk->isotope_frac); - Nread += sscanf(inputLine+117, "%lf", &rlk->hyperfine_frac); - rlk->hyperfine_frac = POW10(rlk->hyperfine_frac); - Nread += sscanf(inputLine+123, "%5d%5d", &hfs_i, &hfs_j); - rlk->hfs_i = ((double) hfs_i) * MILLI * KBOLTZMANN; - rlk->hfs_j = ((double) hfs_j) * MILLI * KBOLTZMANN; - - /* --- Effective Lande factors -- -------------- */ - - Nread += sscanf(inputLine+143, "%5d%5d", &gL_i, &gL_j); - rlk->gL_i = gL_i * MILLI; - rlk->gL_j = gL_j * MILLI; - if (swap_levels) { - SWAPDOUBLE(rlk->hfs_i, rlk->hfs_j); - SWAPDOUBLE(rlk->gL_i, rlk->gL_j); - } - - /* Nread += sscanf(inputLine+154, "%d", &iso_dl); */ - iso_dl = 0; - rlk->iso_dl = iso_dl * MILLI * ANGSTROM_TO_NM; - - checkNread(Nread, Nrequired=17, routineName, checkPoint=1); - /* - printf(" Line: %f (vacuum), %f (air)\n" - " gi, gj: %f, %f\n" - " Ei, Ej: %e, %e\n" - " Aji: %e\n" - " Grad, GStark, GvdWaals: %e, %e, %e\n" - " VdWaals: %d\n" - " hyperfine_frac, isotope_frac: %f, %f\n" - " cross, alpha: %e, %e\n\n", - rlk->lambda0, lambda_air, - rlk->gi, rlk->gj, rlk->Ei, rlk->Ej, rlk->Aji, - rlk->Grad, rlk->GStark, rlk->GvdWaals, - rlk->vdwaals, - rlk->hyperfine_frac, rlk->isotope_frac, - rlk->cross, rlk->alpha); - */ - rlk++; + getUnsoldcross(rlk); + } + /* --- Radiative broadening -- -------------- */ + + if (Grad != 0.0) { + rlk->Grad = POW10(Grad); + } else { + + /* --- Just take the Einstein Aji value-- -------------- */ + + rlk->Grad = rlk->Aji; + } + /* --- Isotope and hyperfine fractions and slpittings -- ---- */ + + Nread += sscanf(inputLine+106, "%d", &rlk->isotope); + Nread += sscanf(inputLine+108, "%lf", &rlk->isotope_frac); + rlk->isotope_frac = POW10(rlk->isotope_frac); + Nread += sscanf(inputLine+117, "%lf", &rlk->hyperfine_frac); + rlk->hyperfine_frac = POW10(rlk->hyperfine_frac); + Nread += sscanf(inputLine+123, "%5d%5d", &hfs_i, &hfs_j); + rlk->level_i.hfs = ((double) hfs_i) * MILLI * KBOLTZMANN; + rlk->level_j.hfs = ((double) hfs_j) * MILLI * KBOLTZMANN; + + /* --- Effective Lande factors -- -------------- */ + + Nread += sscanf(inputLine+143, "%5d%5d", &gL_i, &gL_j); + rlk->level_i.gL = gL_i * MILLI; + rlk->level_j.gL = gL_j * MILLI; + if (swap_levels) { + SWAPDOUBLE(rlk->level_i.hfs, rlk->level_j.hfs); + SWAPDOUBLE(rlk->level_i.gL, rlk->level_j.gL); + } + + /* Nread += sscanf(inputLine+154, "%d", &iso_dl); */ + iso_dl = 0; + rlk->iso_dl = iso_dl * MILLI * ANGSTROM_TO_NM; + + checkNread(Nread, Nrequired=17, routineName, checkPoint=1); + rlk++; } } + sprintf(messageStr, "Read %d Kurucz lines from file %s\n", Nline, listName); Error(MESSAGE, routineName, messageStr); @@ -454,7 +446,7 @@ flags rlk_opacity(double lambda, int nspect, int mu, bool_t to_obs, bool_t contributes, hunt; int Nwhite, Nblue, Nred, NrecStokes; - double dlamb_wing, *pf, dlamb_char, hc_la, ni_gi, nj_gj, lambda0, kT, + double dlamb_wing, *pf, dlamb_char, hc_la, ni_gi, nj_gj, kT, Bijhc_4PI, twohnu3_c2, hc, fourPI, hc_4PI, *eta_Q, *eta_U, *eta_V, eta_l, *chi_Q, *chi_U, *chi_V, chi_l, *chip_Q, *chip_U, *chip_V, @@ -549,103 +541,104 @@ flags rlk_opacity(double lambda, int nspect, int mu, bool_t to_obs, stage, and if abundance is set -- -------------- */ if ((rlk->stage < element->Nstage - 1) && element->abundance_set) { - contributes = TRUE; - if ((metal = element->model) != NULL) { + contributes = TRUE; + if ((metal = element->model) != NULL) { /* --- If an explicit atomic model is present check that we do not already account for this line in this way - - */ - for (kr = 0; kr < metal->Nline; kr++) { - line = metal->line + kr; - dlamb_wing = line->lambda0 * line->qwing * - (atmos.vmicro_char / CLIGHT); - if (fabs(lambda - line->lambda0) <= dlamb_wing && - metal->stage[line->i] == rlk->stage) { - contributes = FALSE; - break; - } - } - } + for (kr = 0; kr < metal->Nline; kr++) { + line = metal->line + kr; + dlamb_wing = line->lambda0 * line->qwing * + (atmos.vmicro_char / CLIGHT); + if (fabs(lambda - line->lambda0) <= dlamb_wing && + metal->stage[line->i] == rlk->stage) { + contributes = FALSE; + break; + } + } + } } else - contributes = FALSE; + contributes = FALSE; /* --- Get opacity from line -- -------------- */ if (contributes) { - hc_la = (HPLANCK * CLIGHT) / (rlk->lambda0 * NM_TO_M); - Bijhc_4PI = hc_4PI * rlk->Bij * rlk->isotope_frac * - rlk->hyperfine_frac * rlk->gi; - twohnu3_c2 = rlk->Aji / rlk->Bji; - - if (input.rlkscatter) { - if (rlk->stage == 0) { - x = 0.68; - C3 = C / (C2_atom * SQ(rlk->lambda0 * NM_TO_M)); - } else { - x = 0.0; - C3 = C / (C2_ion * SQ(rlk->lambda0 * NM_TO_M)); - } - - dE = rlk->Ej - rlk->Ei; - } + hc_la = (HPLANCK * CLIGHT) / (rlk->lambda0 * NM_TO_M); + Bijhc_4PI = hc_4PI * rlk->Bij * rlk->isotope_frac * + rlk->hyperfine_frac * rlk->level_i.g; + twohnu3_c2 = rlk->Aji / rlk->Bji; + + if (input.rlkscatter) { + if (rlk->stage == 0) { + x = 0.68; + C3 = C / (C2_atom * SQ(rlk->lambda0 * NM_TO_M)); + } else { + x = 0.0; + C3 = C / (C2_ion * SQ(rlk->lambda0 * NM_TO_M)); + } + + dE = rlk->level_j.E - rlk->level_i.E; + } /* --- Set flag that line is present at this wavelength -- -- */ - backgrflags.hasline = TRUE; - if (rlk->polarizable) { - backgrflags.ispolarized = TRUE; - if (rlk->zm == NULL) rlk->zm = RLKZeeman(rlk); - } + backgrflags.hasline = TRUE; + if (rlk->polarizable) { + backgrflags.ispolarized = TRUE; + if (rlk->zm == NULL) RLKZeeman(rlk); + } if (element->n == NULL) { - element->n = matrix_double(element->Nstage, atmos.Nspace); - LTEpops_elem(element); - } + element->n = matrix_double(element->Nstage, atmos.Nspace); + LTEpops_elem(element); + } Linear(atmos.Npf, atmos.Tpf, element->pf[rlk->stage], - atmos.Nspace, atmos.T, pf, hunt=TRUE); + atmos.Nspace, atmos.T, pf, hunt=TRUE); - for (k = 0; k < atmos.Nspace; k++) { - phi = RLKProfile(rlk, k, mu, to_obs, lambda, - &phi_Q, &phi_U, &phi_V, - &psi_Q, &psi_U, &psi_V); + for (k = 0; k < atmos.Nspace; k++) { + phi = RLKProfile(rlk, k, mu, to_obs, lambda, + &phi_Q, &phi_U, &phi_V, + &psi_Q, &psi_U, &psi_V); - if (phi){ - kT = 1.0 / (KBOLTZMANN * atmos.T[k]); - ni_gi = element->n[rlk->stage][k] * exp(-rlk->Ei*kT - pf[k]); + if (phi){ + kT = 1.0 / (KBOLTZMANN * atmos.T[k]); + ni_gi = element->n[rlk->stage][k] * + exp(-rlk->level_i.E * kT - pf[k]); nj_gj = ni_gi * exp(-hc_la * kT); - chi_l = Bijhc_4PI * (ni_gi - nj_gj); - eta_l = Bijhc_4PI * twohnu3_c2 * nj_gj; + chi_l = Bijhc_4PI * (ni_gi - nj_gj); + eta_l = Bijhc_4PI * twohnu3_c2 * nj_gj; - if (input.rlkscatter) { - epsilon = 1.0 / (1.0 + C3 * pow(atmos.T[k], 1.5) / - (atmos.ne[k] * - pow(KBOLTZMANN * atmos.T[k] / dE, 1 + x))); + if (input.rlkscatter) { + epsilon = 1.0 / (1.0 + C3 * pow(atmos.T[k], 1.5) / + (atmos.ne[k] * + pow(KBOLTZMANN * atmos.T[k] / dE, 1 + x))); scatt[k] += (1.0 - epsilon) * chi_l * phi; - chi_l *= epsilon; + chi_l *= epsilon; eta_l *= epsilon; - } - - chi[k] += chi_l * phi; - eta[k] += eta_l * phi; - - if (rlk->zm != NULL && rlk->Grad) { - chi_Q[k] += chi_l * phi_Q; - chi_U[k] += chi_l * phi_U; - chi_V[k] += chi_l * phi_V; - - eta_Q[k] += eta_l * phi_Q; - eta_U[k] += eta_l * phi_U; - eta_V[k] += eta_l * phi_V; - - if (input.magneto_optical) { - chip_Q[k] += chi_l * psi_Q; - chip_U[k] += chi_l * psi_U; - chip_V[k] += chi_l * psi_V; + } + + chi[k] += chi_l * phi; + eta[k] += eta_l * phi; + + if (rlk->zm != NULL && rlk->Grad) { + chi_Q[k] += chi_l * phi_Q; + chi_U[k] += chi_l * phi_U; + chi_V[k] += chi_l * phi_V; + + eta_Q[k] += eta_l * phi_Q; + eta_U[k] += eta_l * phi_U; + eta_V[k] += eta_l * phi_V; + + if (input.magneto_optical) { + chip_Q[k] += chi_l * psi_Q; + chip_U[k] += chi_l * psi_U; + chip_V[k] += chi_l * psi_V; + } + } + } } - } - } - } } } } @@ -721,16 +714,16 @@ double RLKProfile(RLK_Line *rlk, int k, int mu, bool_t to_obs, switch (rlk->zm->q[nz]) { case -1: - phi_sm += rlk->zm->strength[nz] * H; - psi_sm += rlk->zm->strength[nz] * F; - break; + phi_sm += rlk->zm->strength[nz] * H; + psi_sm += rlk->zm->strength[nz] * F; + break; case 0: - phi_pi += rlk->zm->strength[nz] * H; - psi_pi += rlk->zm->strength[nz] * F; - break; + phi_pi += rlk->zm->strength[nz] * H; + psi_pi += rlk->zm->strength[nz] * F; + break; case 1: - phi_sp += rlk->zm->strength[nz] * H; - psi_sp += rlk->zm->strength[nz] * F; + phi_sp += rlk->zm->strength[nz] * H; + psi_sp += rlk->zm->strength[nz] * F; } } phi_sigma = phi_sp + phi_sm; @@ -759,7 +752,7 @@ double RLKProfile(RLK_Line *rlk, int k, int mu, bool_t to_obs, /* ------- begin -------------------------- RLKZeeman.c ------------- */ -ZeemanMultiplet* RLKZeeman(RLK_Line *rlk) +void RLKZeeman(RLK_Line *rlk) { const char routineName[] = "RLKZeeman"; @@ -785,9 +778,12 @@ ZeemanMultiplet* RLKZeeman(RLK_Line *rlk) -- -------------- */ - Jl = (rlk->gi - 1.0) / 2.0; - Ju = (rlk->gj - 1.0) / 2.0; - zm = (ZeemanMultiplet *) malloc(sizeof(ZeemanMultiplet)); + Jl = rlk->level_i.J; + Ju = rlk->level_j.J; + + rlk->zm = (ZeemanMultiplet *) malloc(sizeof(ZeemanMultiplet)); + initZeeman(rlk->zm); + zm = rlk->zm; /* --- Count the number of components -- -------------- */ @@ -805,8 +801,8 @@ ZeemanMultiplet* RLKZeeman(RLK_Line *rlk) /* --- Fill the structure and normalize the strengths -- -------- */ - gLl = Lande(rlk->Si, rlk->Li, Jl); - gLu = Lande(rlk->Sj, rlk->Lj, Ju); + gLl = RLKLande(&rlk->level_i); + gLu = RLKLande(&rlk->level_j); n = 0; for (Ml = -Jl; Ml <= Jl; Ml++) { @@ -817,7 +813,7 @@ ZeemanMultiplet* RLKZeeman(RLK_Line *rlk) zm->strength[n] = ZeemanStrength(Ju, Mu, Jl, Ml); norm[zm->q[n]+1] += zm->strength[n]; - if (zm->q[n] == 1) g_eff += zm->shift[n] * zm->strength[n]; + if (zm->q[n] == 1) g_eff += zm->shift[n] * zm->strength[n]; n++; } } @@ -825,59 +821,190 @@ ZeemanMultiplet* RLKZeeman(RLK_Line *rlk) for (n = 0; n < zm->Ncomponent; n++) zm->strength[n] /= norm[zm->q[n]+1]; g_eff /= norm[2]; - - return zm; + zm->g_eff = g_eff; } /* ------- end ---------------------------- RLKZeeman.c ------------- */ -/* ------- begin -------------------------- RLKdeterminate.c -------- */ +/* ------- begin -------------------------- RLKdet_level.c ---------- */ + +#define LS_COUNT 2 +#define JK_COUNT 6 +#define JJ_COUNT 4 -bool_t RLKdeterminate(char *labeli, char *labelj, RLK_Line *rlk) +bool_t RLKdet_level(char* label, RLK_level *level) { - const char routineName[] = "RLKZeeman"; + const char routineName[] = "RLKdet_level"; - char **words, orbit[2]; - bool_t invalid; - int count, multiplicity, length, Nread, Ji, Jj; + char **words, orbit[2], Lchar, L1char, lchar, Kchar, l1char, l2char; + char delims_i[] = "({["; + char delims_f[] = ")}]"; + char *ptr_i, *ptr_f, *quanta_str; + bool_t counterror; + int count, multiplicity, length, Nread, M1; - /* --- Get spin and orbital quantum numbers from level labels -- -- */ + /* --- Get spin and orbital quantum numbers from level labels -- - words = getWords(labeli, " ", &count); - if (words[0]) { - length = strlen(words[count-1]); - Nread = sscanf(words[count-1] + length-2, "%d%1s", - &multiplicity, orbit); - free(words); - if (Nread != 2 || !isupper(orbit[0])) return FALSE; + For explicit LS_COUPLING, or JK_ and JJ_COUPLING provide labels with + explicit quantum numbers as follows (note the delimiters, no spaces). + In the Kurucz line list file (note J is always explicitly listed). - rlk->Li = getOrbital(orbit[0]); - rlk->Si = (multiplicity - 1) / 2.0; - Ji = (rlk->gi - 1.0) / 2.0; - } else - return FALSE; + In keyword.input set RLK_EXPLICIT = TRUE + You CAN NOT mix explicit and non-explicit label modes! - words = getWords(labelj, " ", &count); - if (words[0]) { - length = strlen(words[count-1]); - Nread = sscanf(words[count-1] + length-2, "%d%1s", - &multiplicity, orbit); - free(words); - if (Nread != 2 || !isupper(orbit[0])) return FALSE; + See: Landi degl'Innocenti & Landolfi 2004, pp 76-77 + Note: the label HAS to be 10 characters long, and no other + elements of the line in the .kur should be moved!! - rlk->Lj = getOrbital(orbit[0]); - rlk->Sj = (multiplicity - 1) / 2.0; - Jj = (rlk->gj - 1.0) / 2.0; - } else - return FALSE; + LS_COUPLING: [m,L] Exmpl: '[2P] ' --> S = 0.5, L = 1 + JK_COUPLING: (M1L1J1)lmK Exmpl: '(6D4.5)f2K' --> S1 = 2.5, L1 = 2, + J1 = 4.5, l = 3, K = 3.5 + JJ_COUPLING: {j1l1j2l2} Exmpl: '{1.5p2.5s}' --> j1 = 1.5, + l1 = 1, j2 = 2.5, l2 = 0 + + Example for the FeI 1565.2874 line (JK_COUPLING in upper level): + +1565.2874 -0.476 26.00 50377.905 5.0 s6D)4d f7D 56764.763 4.0 s6D9/4f[3] 8.44 -5.00 -7.70K94 0 0 0 0.000 0 0.000 0 0 1510 1542 + + becomes: + +1565.2874 -0.476 26.00 50377.905 5.0 [7D] 56764.763 4.0 (6D4.5)f2k 8.44 -5.00 -7.70K94 0 0 0 0.000 0 0.000 0 0 1510 1542 + + -- -------------- */ + + if (!input.RLK_explicit) { + + words = getWords(label, " ", &count); + if (words[0]) { + length = strlen(words[count-1]); + Nread = sscanf(words[count-1] + length-2, "%d%1s", + &multiplicity, orbit); + free(words); + if (Nread != 2 || !isupper(orbit[0])) return FALSE; + + level->L = getOrbital(orbit[0]); + level->S = (multiplicity - 1) / 2.0; + + level->cpl = LS_COUPLING; + return TRUE; + } else + return FALSE; + } else { + + /* --- Check for presence of any of three allowed delimiters -- - */ + + ptr_i = strpbrk(label, delims_i); + ptr_f = strpbrk(label, delims_f); + + if (ptr_i == NULL || ptr_f == NULL) { + if (ptr_i != NULL && ptr_f == NULL) { + sprintf(messageStr, " Malformed label: missing ending " + "delimiter: %c, label: %s\n", + delims_f[strchr(delims_i, ptr_i[0]) - delims_i], label); + Error(ERROR_LEVEL_2, routineName, messageStr); + } + if (ptr_i == NULL && ptr_f != NULL) { + sprintf(messageStr, " Malformed label: missing beginning " + "delimiter: %c, label: %s\n", + delims_i[strchr(delims_f, ptr_f[0]) - delims_f], label); + Error(ERROR_LEVEL_2, routineName, messageStr); + } + return FALSE; + } + if (ptr_i[0] != delims_i[strchr(delims_f, ptr_f[0]) - delims_f]) { + + /* --- When delimiters are not matching -- -------------- */ + + sprintf(messageStr, + " Malformed label: mismatched delimiters: '%c %c', %s\n", + ptr_i[0], ptr_f[0], label); + Error(ERROR_LEVEL_2, routineName, messageStr); + return FALSE; + } + + length = ptr_f - ptr_i + 2; + quanta_str = (char *) malloc(length); + strncpy(quanta_str, ptr_i, length-1); + quanta_str[length-1] = '\0'; + + counterror = FALSE; + + switch (ptr_i[0]) { + case '[': + if (Nread = sscanf(quanta_str, "[%1d%1c]", + &multiplicity, &Lchar) != LS_COUNT) + counterror = TRUE; + else { + level->S = (multiplicity - 1) / 2.0; + level->L = getOrbital(Lchar); + + level->cpl = LS_COUPLING; + } + break; - /* --- For the moment only allow electronic dipole transitions -- --*/ + case '(': + if (Nread = sscanf(label, "(%1d%1c%3lf)%1c%1d%1c", + &M1, &L1char, &level->J1, &lchar, + &multiplicity, &Kchar) != JK_COUNT) + counterror = TRUE; + else { + level->S1 = (M1 - 1) / 2.0; + level->L1 = getOrbital(L1char); + level->l = getOrbital(toupper(lchar)); + level->K = getJK_K(Kchar); + + level->cpl = JK_COUPLING; + } + break; - /* if (fabs(Ji - Jj) > 1.0) - return FALSE; - else */ + case '{': + if (Nread = sscanf(quanta_str, "{%3lf%1c%3lf%1c}", + &level->j1, &l1char, + &level->j2, &l2char) != JJ_COUNT) + counterror = TRUE; + else { + level->l1 = getOrbital(toupper(l1char)); + level->l2 = getOrbital(toupper(l2char)); + + level->cpl = JJ_COUPLING; + } + break; + } + if (counterror) { + sprintf(messageStr, "Wrong quantum number count: %s: %d\n", + quanta_str, Nread); + Error(ERROR_LEVEL_2, routineName, messageStr); + return FALSE; + } + free(quanta_str); return TRUE; + } +} +/* ------- end ---------------------------- RLKdet_level.c ---------- */ + +/* ------- begin -------------------------- RLKdet_level_ac.c ---------- */ + +bool_t RLKdet_level_ac(char* label, RLK_level *level) +{ + const char routineName[] = "RLKdet_level_ac"; + + char **words, orbit[2]; + int count, length, Nread; + + /* --- Get orbital quantum numbers from atomic configuration, + adapted from JdlCR STiC -- */ + + /* --- Get spin and orbital quantum numbers from level labels -- -- */ + words = getWords(label, " ", &count); + if (words[0]) { + length = strlen(words[0]); + Nread = sscanf(words[0] + length-1, "%1s", orbit); + free(words); + level->L_ac = getOrbital(toupper(orbit[0])); + } else return FALSE; + + return TRUE; } -/* ------- end ---------------------------- RLKdeterminate.c -------- */ +/* ------- end ---------------------------- RLKdet_level_ac.c ---------- */ /* ------- begin -------------------------- initRLK.c --------------- */ @@ -906,8 +1033,8 @@ void getUnsoldcross(RLK_Line *rlk) } Z = rlk->stage + 1; - deltaR = SQ(E_RYDBERG/(element->ionpot[rlk->stage] - rlk->Ej)) - - SQ(E_RYDBERG/(element->ionpot[rlk->stage] - rlk->Ei)); + deltaR = SQ(E_RYDBERG/(element->ionpot[rlk->stage] - rlk->level_j.E)) - + SQ(E_RYDBERG/(element->ionpot[rlk->stage] - rlk->level_i.E)); if (deltaR <= 0.0) { rlk->vdwaals = KURUCZ; @@ -939,3 +1066,72 @@ void free_BS(Barklemstruct *bs) freeMatrix((void **) bs->alpha); } /* ------- end ---------------------------- free_BS.c --------------- */ + +/* ------- begin -------------------------- RLKLande.c -------------- */ + +double RLKLande(RLK_level *level) +{ + const char routineName[] = "RLKLande"; + + /* --- Lande g factors for different angular momentum coupling + schemes. + + See: Landi degl'Innocenti & Landolfi 2004, pp 76-77 -- ----- */ + + switch (level->cpl) { + case LS_COUPLING: + level->gL = 1.0 + zm_gamma(level->J, level->S, level->L); + + break; + + case JK_COUPLING: + level->gL = 1.0 + zm_gamma(level->J, 0.5, level->K) + + zm_gamma(level->J, level->K, 0.5) * + zm_gamma(level->K, level->J1, level->l) * + zm_gamma(level->J1, level->S1, level->L1); + + break; + + case JJ_COUPLING: + level->gL = 1.0 + zm_gamma(level->J, level->j1, level->j2) * + zm_gamma(level->j1, 0.5, level->l1) + + zm_gamma(level->J, level->j2, level->j1) * + zm_gamma(level->j2, 0.5, level->l2); + + break; + + default: + sprintf(messageStr, "Invalid coupling: %d", level->cpl); + Error(ERROR_LEVEL_2, routineName, messageStr); + } + return level->gL; +} +/* ------- end ---------------------------- RLKLande.c -------------- */ + +/* ------- begin -------------------------- getJK_K.c -- ------------ */ + +double getJK_K(char Kchar) +{ + const char routineName[] = "getJK_K"; + + double K = 0.0; + + switch (Kchar) { + case 'p': K = 0.5; break; + case 'f': K = 1.5; break; + case 'h': K = 2.5; break; + case 'k': K = 3.5; break; + case 'm': K = 4.5; break; + case 'o': K = 5.5; break; + case 'r': K = 6.5; break; + case 't': K = 7.5; break; + case 'u': K = 8.5; break; + case 'v': K = 9.5; break; + case 'w': K = 10.5; break; + default: + sprintf(messageStr, "Invalid Kchar: %c", Kchar); + Error(ERROR_LEVEL_2, routineName, messageStr); + } + return K; +} +/* ------- end ---------------------------- getJK_K.c --------------- */ diff --git a/ltepops.c b/ltepops.c index e8200e0..1191006 100644 --- a/ltepops.c +++ b/ltepops.c @@ -233,8 +233,6 @@ void SetLTEQuantities(void) /* --- Read the collisional data (in MULTI's GENCOL format). After this we can close the input file for the active atom. -- -------------- */ - // Tiago: this must work with atom as string, must use - // atom->offset_coll, and somehow read the file again! (filename gone) CollisionRate(atom, atom->offset_coll); /* --- Compute the fixed rates and store in Cij -- ------------ */ diff --git a/molzeeman.c b/molzeeman.c index 9392b2e..339466f 100644 --- a/molzeeman.c +++ b/molzeeman.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Thu Jun 30 15:35:10 2011 -- + Last modified: Fri May 14 19:15:05 2021 -- -------------------------- ----------RH-- */ @@ -38,7 +38,7 @@ double MolZeemanStr(double Ju, double Mu, double Jl, double Ml) const char routineName[] = "MolZeemanStr"; int q, dJ; - double s; + double s = 0.0; /* --- Return the strength of Zeeman component (Ju, Mu) -> (Jl, Ml), where J and M are the total angular momentum and magnetic @@ -195,7 +195,7 @@ double MolLande_eff(MolecularLine *mrt) /* ------- begin -------------------------- MolZeeman.c ------------- */ -ZeemanMultiplet* MolZeeman(MolecularLine *mrt) +void MolZeeman(MolecularLine *mrt) { const char routineName[] = "MolZeeman"; @@ -221,7 +221,9 @@ ZeemanMultiplet* MolZeeman(MolecularLine *mrt) -- -------------- */ - zm = (ZeemanMultiplet *) malloc(sizeof(ZeemanMultiplet)); + mrt->zm = (ZeemanMultiplet *) malloc(sizeof(ZeemanMultiplet)); + initZeeman(mrt->zm); + zm = mrt->zm; if (mrt->g_Lande_eff != 0.0) { @@ -315,7 +317,5 @@ ZeemanMultiplet* MolZeeman(MolecularLine *mrt) "Zeeman components, gL_eff = %7.4f\n", mrt->molecule->ID, lambda_air, zm->Ncomponent, mrt->g_Lande_eff); Error(MESSAGE, routineName, messageStr); - - return zm; } -/* ------- end ---------------------------- Zeeman.c ---------------- */ +/* ------- end ------------------------- MolZeeman.c ---------------- */ diff --git a/opacity.c b/opacity.c index 6d6e04d..0806013 100644 --- a/opacity.c +++ b/opacity.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Fri Jul 24 16:39:47 2009 -- + Last modified: Fri May 14 17:53:23 2021 -- -------------------------- ----------RH-- */ @@ -444,9 +444,9 @@ void Opacity(int nspect, int mu, bool_t to_obs, bool_t initialize) void alloc_as(int nspect, bool_t crosscoupling) { - register int n, m, nact; + register int m, nact; - int i, j, NrecStokes, NrecStokes_as, Nactive, nt; + int i, j, NrecStokes, NrecStokes_as, nt; Atom *atom; Molecule *molecule; ActiveSet *as; @@ -547,7 +547,7 @@ void alloc_as(int nspect, bool_t crosscoupling) void free_as(int nspect, bool_t crosscoupling) { - register int nact, n, m; + register int nact, m; int i, j, nt; Atom *atom; @@ -582,23 +582,23 @@ void free_as(int nspect, bool_t crosscoupling) freeMatrix((void **) atom->rhth[nt].wla); if (crosscoupling) { - for (m = 0; m < as->Nlower[nact]; m++) { - i = as->lower_levels[nact][m]; - free(atom->rhth[nt].chi_up[i]); - atom->rhth[nt].chi_up[i] = NULL; - } - free(atom->rhth[nt].chi_up); - - for (m = 0; m < as->Nupper[nact]; m++) { - j = as->upper_levels[nact][m]; - free(atom->rhth[nt].chi_down[j]); - free(atom->rhth[nt].Uji_down[j]); - - atom->rhth[nt].chi_down[j] = NULL; - atom->rhth[nt].Uji_down[j] = NULL; - } - free(atom->rhth[nt].chi_down); - free(atom->rhth[nt].Uji_down); + for (m = 0; m < as->Nlower[nact]; m++) { + i = as->lower_levels[nact][m]; + free(atom->rhth[nt].chi_up[i]); + atom->rhth[nt].chi_up[i] = NULL; + } + free(atom->rhth[nt].chi_up); + + for (m = 0; m < as->Nupper[nact]; m++) { + j = as->upper_levels[nact][m]; + free(atom->rhth[nt].chi_down[j]); + free(atom->rhth[nt].Uji_down[j]); + + atom->rhth[nt].chi_down[j] = NULL; + atom->rhth[nt].Uji_down[j] = NULL; + } + free(atom->rhth[nt].chi_down); + free(atom->rhth[nt].Uji_down); } } } @@ -663,16 +663,16 @@ bool_t containsBoundBound(ActiveSet *as) bool_t containsActive(ActiveSet *as) { - register int n, nact; + register int nact; for (nact = 0; nact < atmos.Nactiveatom; nact++) { if (as->Nactiveatomrt[nact] > 0) - return TRUE; + return TRUE; } for (nact = 0; nact < atmos.Nactivemol; nact++) { if (as->Nactivemolrt[nact] > 0) - return TRUE; + return TRUE; } return FALSE; @@ -688,8 +688,8 @@ bool_t containsPRDline(ActiveSet *as) for (nact = 0; nact < atmos.Nactiveatom; nact++) { for (n = 0; n < as->Nactiveatomrt[nact]; n++) { if (as->art[nact][n].type == ATOMIC_LINE && - as->art[nact][n].ptype.line->PRD) { - return TRUE; + as->art[nact][n].ptype.line->PRD) { + return TRUE; } } } @@ -837,7 +837,7 @@ flags MolecularOpacity(double lambda, int nspect, int mu, bool_t to_obs, backgrflags.hasline = TRUE; if (mrt->polarizable) { backgrflags.ispolarized = TRUE; - if (mrt->zm == NULL) mrt->zm = MolZeeman(mrt); + if (mrt->zm == NULL) MolZeeman(mrt); } for (k = 0; k < atmos.Nspace; k++) { diff --git a/profile.c b/profile.c index 9323734..7afc8e4 100644 --- a/profile.c +++ b/profile.c @@ -52,8 +52,6 @@ /* --- Function prototypes -- -------------- */ -void freeZeeman(ZeemanMultiplet *zm); - /* --- Global variables -- -------------- */ @@ -78,14 +76,6 @@ void Profile(AtomicLine *line) *phi_V = NULL, *psi_Q = NULL, *psi_U = NULL, *psi_V = NULL; Atom *atom = line->atom; - ZeemanMultiplet *zm = NULL; - - if (!line->Voigt) { - sprintf(messageStr, - "Magnetic lines cannot have GAUSSian profiles. Line %d -> %d", - line->j, line->i); - Error(ERROR_LEVEL_2, routineName, messageStr); - } getCPU(3, TIME_START, NULL); @@ -106,13 +96,13 @@ void Profile(AtomicLine *line) line->rho_prd = matrix_double(Nlamu, atmos.Nspace); for (la = 0; la < Nlamu; la++) { for (k = 0; k < atmos.Nspace; k++) - line->rho_prd[la][k] = 1.0; + line->rho_prd[la][k] = 1.0; } line->Qelast = (double *) malloc(atmos.Nspace * sizeof(double)); } vbroad = atom->vbroad; - adamp = (double *) malloc(atmos.Nspace * sizeof(double)); + adamp = (double *) calloc(atmos.Nspace, sizeof(double)); if (line->Voigt) Damping(line, adamp); line->wphi = (double *) calloc(atmos.Nspace, sizeof(double)); @@ -120,10 +110,10 @@ void Profile(AtomicLine *line) if (line->polarizable && (input.StokesMode > FIELD_FREE)) { Larmor = (Q_ELECTRON / (4.0*PI*M_ELECTRON)) * (line->lambda0*NM_TO_M); - zm = Zeeman(line); + Zeeman(line); sprintf(messageStr, " -- Atom %2s, line %3d -> %3d has %2d Zeeman components\n", - atom->ID, line->j, line->i, zm->Ncomponent); + atom->ID, line->j, line->i, line->zm->Ncomponent); Error(MESSAGE, routineName, messageStr); } @@ -134,12 +124,12 @@ void Profile(AtomicLine *line) "profile.%.1s_%d-%d.dat" : "profile.%.2s_%d-%d.dat", atom->ID, line->j, line->i); if ((line->fd_profile = - open(filename, O_RDWR | O_CREAT, PERMISSIONS)) == -1) { + open(filename, O_RDWR | O_CREAT, PERMISSIONS)) == -1) { sprintf(messageStr, "Unable to open profile file %s", filename); Error(ERROR_LEVEL_2, routineName, messageStr); } - if (line->polarizable && (input.StokesMode > FIELD_FREE)) { + if (line->polarizable && (input.StokesMode == FULL_STOKES)) { NrecStokes = (input.magneto_optical) ? 7 : 4; phi = (double *) malloc(NrecStokes*atmos.Nspace * sizeof(double)); @@ -150,9 +140,9 @@ void Profile(AtomicLine *line) phi_V = phi + 3*atmos.Nspace; if (input.magneto_optical) { - psi_Q = phi + 4*atmos.Nspace; - psi_U = phi + 5*atmos.Nspace; - psi_V = phi + 6*atmos.Nspace; + psi_Q = phi + 4*atmos.Nspace; + psi_U = phi + 5*atmos.Nspace; + psi_V = phi + 6*atmos.Nspace; } } else { NrecStokes = 1; @@ -160,26 +150,26 @@ void Profile(AtomicLine *line) } } else { if (atmos.moving || - (line->polarizable && (input.StokesMode > FIELD_FREE))) { + (line->polarizable && (input.StokesMode == FULL_STOKES))) { Nlamu = 2*atmos.Nrays*line->Nlambda; line->phi = matrix_double(Nlamu, atmos.Nspace); - if (line->polarizable && (input.StokesMode > FIELD_FREE)) { - line->phi_Q = matrix_double(Nlamu, atmos.Nspace); - line->phi_U = matrix_double(Nlamu, atmos.Nspace); - line->phi_V = matrix_double(Nlamu, atmos.Nspace); + if (line->polarizable && (input.StokesMode == FULL_STOKES)) { + line->phi_Q = matrix_double(Nlamu, atmos.Nspace); + line->phi_U = matrix_double(Nlamu, atmos.Nspace); + line->phi_V = matrix_double(Nlamu, atmos.Nspace); - if (input.magneto_optical) { - line->psi_Q = matrix_double(Nlamu, atmos.Nspace); - line->psi_U = matrix_double(Nlamu, atmos.Nspace); - line->psi_V = matrix_double(Nlamu, atmos.Nspace); - } + if (input.magneto_optical) { + line->psi_Q = matrix_double(Nlamu, atmos.Nspace); + line->psi_U = matrix_double(Nlamu, atmos.Nspace); + line->psi_V = matrix_double(Nlamu, atmos.Nspace); + } } } else line->phi = matrix_double(line->Nlambda, atmos.Nspace); } - if (line->polarizable && (input.StokesMode > FIELD_FREE)) { + if (line->polarizable && (input.StokesMode == FULL_STOKES)) { /* --- Temporary storage for inner loop variables, vB is the Zeeman splitting due to the local magnetic field -- ------ */ @@ -196,137 +186,141 @@ void Profile(AtomicLine *line) /* --- Calculate the absorption profile and store for each line -- */ if (atmos.moving || - (line->polarizable && (input.StokesMode > FIELD_FREE))) { + (line->polarizable && (input.StokesMode == FULL_STOKES))) { - v_los = matrix_double(atmos.Nrays, atmos.Nspace); - for (mu = 0; mu < atmos.Nrays; mu++) { - for (k = 0; k < atmos.Nspace; k++) { - v_los[mu][k] = vproject(k, mu) / vbroad[k]; + if (atmos.moving) { + v_los = matrix_double(atmos.Nrays, atmos.Nspace); + for (mu = 0; mu < atmos.Nrays; mu++) { + for (k = 0; k < atmos.Nspace; k++) { + v_los[mu][k] = vproject(k, mu) / vbroad[k]; + } } } v = matrix_double(atmos.Nspace, line->Ncomponent); for (la = 0; la < line->Nlambda; la++) { for (n = 0; n < line->Ncomponent; n++) { - for (k = 0; k < atmos.Nspace; k++) { - v[k][n] = (line->lambda[la] - line->lambda0 - line->c_shift[n]) * - CLIGHT / (vbroad[k] * line->lambda0); - } + for (k = 0; k < atmos.Nspace; k++) { + v[k][n] = (line->lambda[la] - line->lambda0 - line->c_shift[n]) * + CLIGHT / (vbroad[k] * line->lambda0); + } } for (mu = 0; mu < atmos.Nrays; mu++) { - wlamu = getwlambda_line(line, la) * 0.5*atmos.wmu[mu]; - - for (to_obs = 0; to_obs <= 1; to_obs++) { - sign = (to_obs) ? 1.0 : -1.0; - lamu = 2*(atmos.Nrays*la + mu) + to_obs; - - /* --- Assign pointers to the proper phi and psi arrays and - zero the profiles in case of conservative memory - option. In the normal case the call matrix_double - initializes the whole array to zero -- ------------- */ - - if (input.limit_memory) { - for (k = 0; k< NrecStokes*atmos.Nspace; k++) phi[k] = 0.0; - } else { - phi = line->phi[lamu]; - if (line->polarizable && (input.StokesMode > FIELD_FREE)) { - phi_Q = line->phi_Q[lamu]; - phi_U = line->phi_U[lamu]; - phi_V = line->phi_V[lamu]; - - if (input.magneto_optical) { - psi_Q = line->psi_Q[lamu]; - psi_U = line->psi_U[lamu]; - psi_V = line->psi_V[lamu]; - } - } - } - - if (line->polarizable && (input.StokesMode > FIELD_FREE)) { - for (k = 0; k < atmos.Nspace; k++) { - sin2_gamma = 1.0 - SQ(atmos.cos_gamma[mu][k]); - - /* --- For the sign conventions to the phi and psi - contributions depending on the direction along the ray - - See: - -- A. van Ballegooijen: "Radiation in Strong Magnetic - Fields", in Numerical Radiative Transfer, W. Kalkofen - 1987, p. 285 -- -------------- */ - + wlamu = getwlambda_line(line, la) * 0.5*atmos.wmu[mu]; + + for (to_obs = 0; to_obs <= 1; to_obs++) { + sign = (to_obs) ? 1.0 : -1.0; + lamu = 2*(atmos.Nrays*la + mu) + to_obs; + + /* --- Assign pointers to the proper phi and psi arrays and + zero the profiles in case of conservative memory + option. In the normal case the call matrix_double + initializes the whole array to zero -- ------------- */ + + if (input.limit_memory) { + for (k = 0; k< NrecStokes*atmos.Nspace; k++) phi[k] = 0.0; + } else { + phi = line->phi[lamu]; + if (line->polarizable && (input.StokesMode == FULL_STOKES)) { + phi_Q = line->phi_Q[lamu]; + phi_U = line->phi_U[lamu]; + phi_V = line->phi_V[lamu]; + + if (input.magneto_optical) { + psi_Q = line->psi_Q[lamu]; + psi_U = line->psi_U[lamu]; + psi_V = line->psi_V[lamu]; + } + } + } + + if (line->polarizable && (input.StokesMode == FULL_STOKES)) { + for (k = 0; k < atmos.Nspace; k++) { + sin2_gamma = 1.0 - SQ(atmos.cos_gamma[mu][k]); + + /* --- For the sign conventions to the phi and psi + contributions depending on the direction along the ray + + See: + -- A. van Ballegooijen: "Radiation in Strong Magnetic + Fields", in Numerical Radiative Transfer, W. Kalkofen + 1987, p. 285 -- -------------- */ + /* --- Sum over isotopes -- -------------- */ - for (n = 0; n < line->Ncomponent; n++) { - vk = v[k][n] + sign * v_los[mu][k]; - - phi_sm = phi_pi = phi_sp = 0.0; - psi_sm = psi_pi = psi_sp = 0.0; - - /* --- Sum over Zeeman sub-levels -- -------------- */ - - for (nz = 0; nz < zm->Ncomponent; nz++) { - H = Voigt(adamp[k], vk - zm->shift[nz]*vB[k], - &F, HUMLICEK); - - switch (zm->q[nz]) { - case -1: - phi_sm += zm->strength[nz] * H; - psi_sm += zm->strength[nz] * F; - break; - case 0: - phi_pi += zm->strength[nz] * H; - psi_pi += zm->strength[nz] * F; - break; - case 1: - phi_sp += zm->strength[nz] * H; - psi_sp += zm->strength[nz] * F; - } - } - phi_sigma = (phi_sp + phi_sm) * line->c_fraction[n]; - phi_delta = 0.5*phi_pi * line->c_fraction[n] - 0.25*phi_sigma; - - phi[k] += (phi_delta*sin2_gamma + 0.5*phi_sigma) * sv[k]; - phi_Q[k] += sign * - phi_delta * sin2_gamma * atmos.cos_2chi[mu][k] * sv[k]; - phi_U[k] += - phi_delta * sin2_gamma * atmos.sin_2chi[mu][k] * sv[k]; - phi_V[k] += sign * - 0.5*(phi_sp - phi_sm) * atmos.cos_gamma[mu][k] * sv[k]; - - if (input.magneto_optical) { - psi_sigma = (psi_sp + psi_sm) * line->c_fraction[n]; - psi_delta = 0.5*psi_pi * line->c_fraction[n] - - 0.25*psi_sigma; - - psi_Q[k] += sign * - psi_delta * sin2_gamma * atmos.cos_2chi[mu][k] * sv[k]; - psi_U[k] += - psi_delta * sin2_gamma * atmos.sin_2chi[mu][k] * sv[k]; - psi_V[k] += sign * - 0.5 * (psi_sp - psi_sm) * atmos.cos_gamma[mu][k] * sv[k]; - } - } + for (n = 0; n < line->Ncomponent; n++) { + vk = v[k][n]; + if (atmos.moving) vk += sign * v_los[mu][k]; + + phi_sm = phi_pi = phi_sp = 0.0; + psi_sm = psi_pi = psi_sp = 0.0; + + /* --- Sum over Zeeman sub-levels -- -------------- */ + + for (nz = 0; nz < line->zm->Ncomponent; nz++) { + H = Voigt(adamp[k], vk - line->zm->shift[nz]*vB[k], + &F, HUMLICEK); + + switch (line->zm->q[nz]) { + case -1: + phi_sm += line->zm->strength[nz] * H; + psi_sm += line->zm->strength[nz] * F; + break; + case 0: + phi_pi += line->zm->strength[nz] * H; + psi_pi += line->zm->strength[nz] * F; + break; + case 1: + phi_sp += line->zm->strength[nz] * H; + psi_sp += line->zm->strength[nz] * F; + } + } + phi_sigma = (phi_sp + phi_sm) * line->c_fraction[n]; + phi_delta = 0.5*phi_pi * line->c_fraction[n] - 0.25*phi_sigma; + + phi[k] += (phi_delta*sin2_gamma + 0.5*phi_sigma) * sv[k]; + phi_Q[k] += sign * + phi_delta * sin2_gamma * atmos.cos_2chi[mu][k] * sv[k]; + phi_U[k] += + phi_delta * sin2_gamma * atmos.sin_2chi[mu][k] * sv[k]; + phi_V[k] += sign * + 0.5*(phi_sp - phi_sm) * atmos.cos_gamma[mu][k] * sv[k]; + + if (input.magneto_optical) { + psi_sigma = (psi_sp + psi_sm) * line->c_fraction[n]; + psi_delta = 0.5*psi_pi * line->c_fraction[n] - + 0.25*psi_sigma; + + psi_Q[k] += sign * + psi_delta * sin2_gamma * atmos.cos_2chi[mu][k] * sv[k]; + psi_U[k] += + psi_delta * sin2_gamma * atmos.sin_2chi[mu][k] * sv[k]; + psi_V[k] += sign * + 0.5 * (psi_sp - psi_sm) * atmos.cos_gamma[mu][k] * sv[k]; + } + } + /* --- Ensure proper normalization of the profile -- -- */ - line->wphi[k] += wlamu * phi[k]; - } - } else { + line->wphi[k] += wlamu * phi[k]; + } + } else { - /* --- Field-free case -- -------------- */ + /* --- Field-free case -- -------------- */ for (k = 0; k < atmos.Nspace; k++) { - for (n = 0; n < line->Ncomponent; n++) { - vk = v[k][n] + sign * v_los[mu][k]; - - phi[k] += Voigt(adamp[k], vk, NULL, ARMSTRONG) * - line->c_fraction[n] / (SQRTPI * atom->vbroad[k]); + for (n = 0; n < line->Ncomponent; n++) { + vk = v[k][n] + sign * v_los[mu][k]; + + phi[k] += Voigt(adamp[k], vk, NULL, ARMSTRONG) * + line->c_fraction[n] / (SQRTPI * atom->vbroad[k]); + } + line->wphi[k] += phi[k] * wlamu; + } + } + if (input.limit_memory) writeProfile(line, lamu, phi); } - line->wphi[k] += phi[k] * wlamu; - } - } - if (input.limit_memory) writeProfile(line, lamu, phi); - } } } } else { @@ -337,18 +331,18 @@ void Profile(AtomicLine *line) wlamu = getwlambda_line(line, la); if (input.limit_memory) - for (k = 0; k < atmos.Nspace; k++) phi[k] = 0.0; + for (k = 0; k < atmos.Nspace; k++) phi[k] = 0.0; else - phi = line->phi[la]; + phi = line->phi[la]; for (k = 0; k < atmos.Nspace; k++) { - for (n = 0; n < line->Ncomponent; n++) { - vk = (line->lambda[la] - line->lambda0 - line->c_shift[n]) * - CLIGHT / (line->lambda0 * atom->vbroad[k]); - phi[k] += Voigt(adamp[k], vk, NULL, ARMSTRONG) * - line->c_fraction[n] / (SQRTPI * atom->vbroad[k]); - } - line->wphi[k] += phi[k] * wlamu; + for (n = 0; n < line->Ncomponent; n++) { + vk = (line->lambda[la] - line->lambda0 - line->c_shift[n]) * + CLIGHT / (line->lambda0 * atom->vbroad[k]); + phi[k] += Voigt(adamp[k], vk, NULL, ARMSTRONG) * + line->c_fraction[n] / (SQRTPI * atom->vbroad[k]); + } + line->wphi[k] += phi[k] * wlamu; } if (input.limit_memory) writeProfile(line, la, phi); } @@ -362,16 +356,20 @@ void Profile(AtomicLine *line) free(adamp); if (input.limit_memory) free(phi); - if (atmos.moving || (line->polarizable && (input.StokesMode > FIELD_FREE))) { - if (zm) - freeZeeman(zm); - free(zm); + if (atmos.moving) freeMatrix((void **) v_los); + if (atmos.moving || + (line->polarizable && (input.StokesMode == FULL_STOKES))) { + if (line->zm != NULL) { + freeZeeman(line->zm); + free(line->zm); + line->zm = NULL; + } + freeMatrix((void **) v); + } + if (line->polarizable && (input.StokesMode == FULL_STOKES)) { free(vB); free(sv); - freeMatrix((void **) v); - freeMatrix((void **) v_los); - sprintf(messageStr, "Stokes prof %7.1f", line->lambda0); } else sprintf(messageStr, "Profile %7.1f", line->lambda0); @@ -424,24 +422,24 @@ void MolecularProfile(MolecularLine *mrt) if (atmos.moving) { for (la = 0; la < mrt->Nlambda; la++) { for (k = 0; k < atmos.Nspace; k++) - v[k] = (mrt->lambda[la] - mrt->lambda0) * CLIGHT / - (molecule->vbroad[k] * mrt->lambda0); + v[k] = (mrt->lambda[la] - mrt->lambda0) * CLIGHT / + (molecule->vbroad[k] * mrt->lambda0); for (mu = 0; mu < atmos.Nrays; mu++) { - wlamu = getwlambda_mrt(mrt, la) * 0.5*atmos.wmu[mu]; + wlamu = getwlambda_mrt(mrt, la) * 0.5*atmos.wmu[mu]; - /* --- First the downward, then the upward direction -- ----- */ + /* --- First the downward, then the upward direction -- ----- */ - phi_down = mrt->phi[2*(atmos.Nrays*la + mu)]; - phi_up = mrt->phi[2*(atmos.Nrays*la + mu) + 1]; + phi_down = mrt->phi[2*(atmos.Nrays*la + mu)]; + phi_up = mrt->phi[2*(atmos.Nrays*la + mu) + 1]; - for (k = 0; k < atmos.Nspace; k++) { - phi_down[k] = Voigt(adamp[k], v[k] - v_los[mu][k], - NULL, ARMSTRONG) * sv[k]; - phi_up[k] = Voigt(adamp[k], v[k] + v_los[mu][k], - NULL, ARMSTRONG) * sv[k]; - mrt->wphi[k] += (phi_down[k] + phi_up[k]) * wlamu; - } + for (k = 0; k < atmos.Nspace; k++) { + phi_down[k] = Voigt(adamp[k], v[k] - v_los[mu][k], + NULL, ARMSTRONG) * sv[k]; + phi_up[k] = Voigt(adamp[k], v[k] + v_los[mu][k], + NULL, ARMSTRONG) * sv[k]; + mrt->wphi[k] += (phi_down[k] + phi_up[k]) * wlamu; + } } } } else { @@ -452,11 +450,11 @@ void MolecularProfile(MolecularLine *mrt) wlambda = getwlambda_mrt(mrt, la); for (k = 0; k < atmos.Nspace; k++) { - vk = (mrt->lambda[la] - mrt->lambda0) * CLIGHT / - (mrt->lambda0 * molecule->vbroad[k]); - mrt->phi[la][k] = Voigt(adamp[k], vk, NULL, ARMSTRONG) * sv[k]; - mrt->wphi[k] += mrt->phi[la][k] * wlambda; - } + vk = (mrt->lambda[la] - mrt->lambda0) * CLIGHT / + (mrt->lambda0 * molecule->vbroad[k]); + mrt->phi[la][k] = Voigt(adamp[k], vk, NULL, ARMSTRONG) * sv[k]; + mrt->wphi[k] += mrt->phi[la][k] * wlambda; + } } } for (k = 0; k < atmos.Nspace; k++) mrt->wphi[k] = 1.0 / mrt->wphi[k]; diff --git a/readatom.c b/readatom.c index 461da8d..7134ea0 100644 --- a/readatom.c +++ b/readatom.c @@ -653,6 +653,7 @@ void initAtomicLine(AtomicLine *line) line->id0 = NULL; line->id1 = NULL; line->gII = NULL; + line->zm = NULL; } /* ------- end ---------------------------- initAtomicLine.c -------- */ @@ -732,6 +733,9 @@ void freeAtomicLine(AtomicLine *line) if (line->psi_Q != NULL) freeMatrix((void **) line->psi_Q); if (line->psi_U != NULL) freeMatrix((void **) line->psi_U); if (line->psi_V != NULL) freeMatrix((void **) line->psi_V); + if (atmos.Stokes && input.StokesMode == FULL_STOKES) { + if (line->zm != NULL)freeZeeman(line->zm); + } } if (line->c_shift != NULL) free(line->c_shift); if (line->c_fraction != NULL) free(line->c_fraction); diff --git a/readinput.c b/readinput.c index 2121a7d..94746ba 100644 --- a/readinput.c +++ b/readinput.c @@ -169,6 +169,9 @@ void readInput(char *input_string) &input.magneto_optical, setboolValue}, {"BACKGROUND_POLARIZATION", "FALSE", FALSE, KEYWORD_DEFAULT, &input.backgr_pol, setboolValue}, + {"RLK_EXPLICIT", "FALSE", FALSE, KEYWORD_DEFAULT, + &input.RLK_explicit, setboolValue}, + {"XDR_ENDIAN", "TRUE", FALSE, KEYWORD_OPTIONAL, &input.xdr_endian, setboolValue}, @@ -317,6 +320,9 @@ void readInput(char *input_string) /* --- Stokes for the moment only in 1D plane -- -------------- */ if (strcmp(input.Stokes_input, "none")) { + + /* --- Magnetic field is specified -- -------------- */ + switch (topology) { case ONE_D_PLANE: case TWO_D_PLANE: @@ -354,6 +360,9 @@ void readInput(char *input_string) "Cannot accomodate magnetic fields in this topology"); } } else { + + /* --- No magnetic field input specified -- -------------- */ + if (atmos.B_char != 0.0) { Error(WARNING, routineName, "Ignoring value of keyword B_STRENGTH_CHAR when no " diff --git a/rh15d/background_p.c b/rh15d/background_p.c index acf646d..7fd9224 100644 --- a/rh15d/background_p.c +++ b/rh15d/background_p.c @@ -223,7 +223,7 @@ void close_Background(void) void Background_p(bool_t write_analyze_output, bool_t equilibria_only) { const char routineName[] = "Background_p"; - register int k, nspect, mu, to_obs; + register int n, k, nspect, mu, to_obs; static int ne_iter = 0; bool_t do_fudge, fromscratch; @@ -276,37 +276,35 @@ void Background_p(bool_t write_analyze_output, bool_t equilibria_only) // sanity check if (input.StokesMode != NO_STOKES) - Error(ERROR_LEVEL_2,routineName,"Polarized transfer and in-memory backgr. op. is not implemented."); + Error(ERROR_LEVEL_2,routineName,"Polarized transfer and in-memory backgr. op. is not implemented."); - if (atmos.chi_b != NULL) freeMatrix((void **) atmos.chi_b); - atmos.chi_b = matrix_double(spectrum.Nspect * 2 * atmos.Nrays, atmos.Nspace); + if (atmos.chi_b != NULL) freeMatrix((void **) atmos.chi_b); + atmos.chi_b = matrix_double(spectrum.Nspect * 2 * atmos.Nrays, atmos.Nspace); - if (atmos.eta_b != NULL) freeMatrix((void **) atmos.eta_b); - atmos.eta_b = matrix_double(spectrum.Nspect * 2 * atmos.Nrays, atmos.Nspace); + if (atmos.eta_b != NULL) freeMatrix((void **) atmos.eta_b); + atmos.eta_b = matrix_double(spectrum.Nspect * 2 * atmos.Nrays, atmos.Nspace); - if (atmos.sca_b != NULL) freeMatrix((void **) atmos.sca_b); - atmos.sca_b = matrix_double(spectrum.Nspect * 2 * atmos.Nrays, atmos.Nspace); + if (atmos.sca_b != NULL) freeMatrix((void **) atmos.sca_b); + atmos.sca_b = matrix_double(spectrum.Nspect * 2 * atmos.Nrays, atmos.Nspace); - } else { - - // on file + } else { - /* --- Open background file -- ------------- */ + /* --- Open background file -- ------------- */ - /* get file name, with the MPI rank */ - sprintf(file_background,"%s_p%d%s", input.background_File , mpi.rank, fext); + /* get file name, with the MPI rank */ + sprintf(file_background,"%s_p%d%s", input.background_File , mpi.rank, fext); - if ((atmos.fd_background = - open(file_background, O_RDWR | O_CREAT | O_TRUNC, PERMISSIONS)) == -1) { - sprintf(messageStr, "Unable to open output file %s", - file_background); - Error(ERROR_LEVEL_2, routineName, messageStr); - } + if ((atmos.fd_background = + open(file_background, O_RDWR | O_CREAT | O_TRUNC, PERMISSIONS)) == -1) { + sprintf(messageStr, "Unable to open output file %s", + file_background); + Error(ERROR_LEVEL_2, routineName, messageStr); + } - /* Zero record number, as background file is being overwritten */ - mpi.backgrrecno = 0; + /* Zero record number, as background file is being overwritten */ + mpi.backgrrecno = 0; - } + } /* --- Allocate temporary storage space. The quantities are used @@ -401,14 +399,14 @@ void Background_p(bool_t write_analyze_output, bool_t equilibria_only) if (Hminus_bf(wavelength, chi, eta)) { for (k = 0; k < atmos.Nspace; k++) { - chi_ai[k] += chi[k]; - eta_ai[k] += eta[k]; + chi_ai[k] += chi[k]; + eta_ai[k] += eta[k]; } } if (Hminus_ff(wavelength, chi)) { for (k = 0; k < atmos.Nspace; k++) { - chi_ai[k] += chi[k]; - eta_ai[k] += chi[k] * Bnu[k]; + chi_ai[k] += chi[k]; + eta_ai[k] += chi[k] * Bnu[k]; } } /* --- Opacity fudge factors, applied to Hminus opacity -- ------ */ @@ -417,30 +415,30 @@ void Background_p(bool_t write_analyze_output, bool_t equilibria_only) Linear(Nfudge, lambda_fudge, fudge[0], 1, &wavelength, &Hmin_fudge, FALSE); for (k = 0; k < atmos.Nspace; k++) { - chi_ai[k] *= Hmin_fudge; - eta_ai[k] *= Hmin_fudge; + chi_ai[k] *= Hmin_fudge; + eta_ai[k] *= Hmin_fudge; } } /* --- Opacities from bound-free transitions in OH and CH -- ---- */ if (OH_bf_opac(wavelength, chi, eta)) { for (k = 0; k < atmos.Nspace; k++) { - chi_ai[k] += chi[k]; - eta_ai[k] += eta[k]; + chi_ai[k] += chi[k]; + eta_ai[k] += eta[k]; } } if (CH_bf_opac(wavelength, chi, eta)) { for (k = 0; k < atmos.Nspace; k++) { - chi_ai[k] += chi[k]; - eta_ai[k] += eta[k]; + chi_ai[k] += chi[k]; + eta_ai[k] += eta[k]; } } /* --- Neutral Hydrogen Bound-Free and Free-Free -- ------------ */ if (Hydrogen_bf(wavelength, chi, eta)) { for (k = 0; k < atmos.Nspace; k++) { - chi_ai[k] += chi[k]; - eta_ai[k] += eta[k]; + chi_ai[k] += chi[k]; + eta_ai[k] += eta[k]; } } Hydrogen_ff(wavelength, chi); @@ -452,22 +450,22 @@ void Background_p(bool_t write_analyze_output, bool_t equilibria_only) if (Rayleigh(wavelength, atmos.H, scatt)) { for (k = 0; k < atmos.Nspace; k++) { - sca_ai[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_ai[k] += scatt[k]; + sca_ai[k] += scatt[k]; } } /* --- Absorption by H + H^+ (referred to as H2plus free-free) -- */ if (H2plus_ff(wavelength, chi)) { for (k = 0; k < atmos.Nspace; k++) { - chi_ai[k] += chi[k]; - eta_ai[k] += chi[k] * Bnu[k]; + chi_ai[k] += chi[k]; + eta_ai[k] += chi[k] * Bnu[k]; } } /* --- Rayleigh scattering and free-free absorption by @@ -475,20 +473,20 @@ void Background_p(bool_t write_analyze_output, bool_t equilibria_only) if (Rayleigh_H2(wavelength, scatt)) { for (k = 0; k < atmos.Nspace; k++) { - sca_ai[k] += scatt[k]; + sca_ai[k] += scatt[k]; } } if (H2minus_ff(wavelength, chi)) { for (k = 0; k < atmos.Nspace; k++) { - chi_ai[k] += chi[k]; - eta_ai[k] += chi[k] * Bnu[k]; + chi_ai[k] += chi[k]; + eta_ai[k] += chi[k] * Bnu[k]; } } /* --- Bound-Free opacities due to ``metals'' -- -------------- */ if (do_fudge) { Linear(Nfudge, lambda_fudge, fudge[2], - 1, &wavelength, &metal_fudge, FALSE); + 1, &wavelength, &metal_fudge, FALSE); } else { metal_fudge = 1.0; } @@ -519,57 +517,55 @@ void Background_p(bool_t write_analyze_output, bool_t equilibria_only) if (atmos.moving || atmos.Stokes) { for (mu = 0; mu < atmos.Nrays; mu++) { for (to_obs = 0; to_obs <= 1; to_obs++) { - index = 2*(nspect*atmos.Nrays + mu) + to_obs; + index = 2*(nspect*atmos.Nrays + mu) + to_obs; /* --- First, copy the angle-independent parts -- --------- */ - 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]; - } + 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) { - for (k = atmos.Nspace; k < 4*atmos.Nspace; k++) { - chi_c[k] = 0.0; - eta_c[k] = 0.0; + if (atmos.Stokes) { + for (k = atmos.Nspace; k < 4*atmos.Nspace; k++) { + chi_c[k] = 0.0; + eta_c[k] = 0.0; } if (input.magneto_optical) for (k = 0; k < 3*atmos.Nspace; k++) chip_c[k] = 0.0; - } + } /* --- Add opacity from passive atomic lines (including hydrogen) -- -------------- */ - if (input.allow_passive_bb) { - backgrflags = passive_bb(wavelength, nspect, mu, to_obs, - chi, eta, chip); - if (backgrflags.hasline) { - atmos.backgrflags[nspect].hasline = TRUE; - if (backgrflags.ispolarized) { + if (input.allow_passive_bb) { + backgrflags = passive_bb(wavelength, nspect, mu, to_obs, + chi, eta, chip); + if (backgrflags.hasline) { + atmos.backgrflags[nspect].hasline = TRUE; + if (backgrflags.ispolarized) { NrecStokes = 4; atmos.backgrflags[nspect].ispolarized = TRUE; - if (input.magneto_optical) { + if (input.magneto_optical) { for (k = 0; k < 3*atmos.Nspace; k++) chip_c[k] += chip[k]; } - } else - NrecStokes = 1; + } else NrecStokes = 1; for (k = 0; k < NrecStokes*atmos.Nspace; k++) { chi_c[k] += chi[k]; eta_c[k] += eta[k]; } - - } - } + } + } /* --- Add opacity from Kurucz line list -- -------------- */ if (atmos.Nrlk > 0) { - backgrflags = rlk_opacity(wavelength, nspect, mu, to_obs, - chi, eta, scatt, chip); - if (backgrflags.hasline) { - atmos.backgrflags[nspect].hasline = TRUE; + backgrflags = rlk_opacity(wavelength, nspect, mu, to_obs, + chi, eta, scatt, chip); + if (backgrflags.hasline) { + atmos.backgrflags[nspect].hasline = TRUE; if (backgrflags.ispolarized) { NrecStokes = 4; atmos.backgrflags[nspect].ispolarized = TRUE; @@ -577,27 +573,26 @@ void Background_p(bool_t write_analyze_output, bool_t equilibria_only) for (k = 0; k < 3*atmos.Nspace; k++) chip_c[k] += chip[k]; } - } else - NrecStokes = 1; + } else NrecStokes = 1; for (k = 0; k < NrecStokes*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 -- -------------- */ - - backgrflags = MolecularOpacity(wavelength, nspect, mu, to_obs, - chi, eta, chip); - if (backgrflags.hasline) { - atmos.backgrflags[nspect].hasline = TRUE; + 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 -- -------------- */ + + backgrflags = MolecularOpacity(wavelength, nspect, mu, to_obs, + chi, eta, chip); + if (backgrflags.hasline) { + atmos.backgrflags[nspect].hasline = TRUE; if (backgrflags.ispolarized) { NrecStokes = 4; atmos.backgrflags[nspect].ispolarized = TRUE; @@ -605,38 +600,37 @@ void Background_p(bool_t write_analyze_output, bool_t equilibria_only) for (k = 0; k < 3*atmos.Nspace; k++) chip_c[k] += chip[k]; } - } else - NrecStokes = 1; + } else NrecStokes = 1; for (k = 0; k < NrecStokes*atmos.Nspace; k++) { chi_c[k] += chi[k]; eta_c[k] += eta[k]; } - } - /* --- Store angle-dependent results only if at least one - line was found at this wavelength -- -------------- */ - - atmos.backgrrecno[index] = mpi.backgrrecno; - if ((mu == atmos.Nrays-1 && to_obs) || - (atmos.backgrflags[nspect].hasline && - (atmos.moving || atmos.backgrflags[nspect].ispolarized))) { - - if (input.backgr_in_mem) { - - if ((mu == atmos.Nrays-1 && to_obs) && !(atmos.backgrflags[nspect].hasline && - (atmos.moving || atmos.backgrflags[nspect].ispolarized)) ){ - storeBackground(nspect, 0, 0, chi_c, eta_c, sca_c); - } else { - storeBackground(nspect, mu, to_obs, chi_c, eta_c, sca_c); - } - - } else { - - mpi.backgrrecno += writeBackground(nspect, mu, to_obs, - chi_c, eta_c, sca_c, chip_c); - } - } - } + } + /* --- Store angle-dependent results only if at least one + line was found at this wavelength -- -------------- */ + + atmos.backgrrecno[index] = mpi.backgrrecno; + if ((mu == atmos.Nrays-1 && to_obs) || + (atmos.backgrflags[nspect].hasline && + (atmos.moving || atmos.backgrflags[nspect].ispolarized))) { + + if (input.backgr_in_mem) { + + if ((mu == atmos.Nrays-1 && to_obs) && !(atmos.backgrflags[nspect].hasline && + (atmos.moving || atmos.backgrflags[nspect].ispolarized)) ){ + storeBackground(nspect, 0, 0, chi_c, eta_c, sca_c); + } else { + storeBackground(nspect, mu, to_obs, chi_c, eta_c, sca_c); + } + + } else { + + mpi.backgrrecno += writeBackground(nspect, mu, to_obs, + chi_c, eta_c, sca_c, chip_c); + } + } + } } } else { @@ -644,55 +638,55 @@ void Background_p(bool_t write_analyze_output, bool_t equilibria_only) atomic lines (including hydrogen) -- -------------- */ if (input.allow_passive_bb) { - backgrflags = passive_bb(wavelength, nspect, 0, TRUE, - chi, eta, 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]; - } - } + backgrflags = passive_bb(wavelength, nspect, 0, TRUE, + chi, eta, 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]; + } + } } /* --- Add opacity from Kurucz line list -- -------------- */ if (atmos.Nrlk > 0) { - backgrflags = rlk_opacity(wavelength, nspect, 0, TRUE, - 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]; - } - } - } + backgrflags = rlk_opacity(wavelength, nspect, 0, TRUE, + 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 -- -------------- */ backgrflags = MolecularOpacity(wavelength, nspect, 0, TRUE, - chi, eta, NULL); + chi, eta, 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]; - } + atmos.backgrflags[nspect].hasline = TRUE; + for (k = 0; k < atmos.Nspace; k++) { + chi_c[k] += chi[k]; + eta_c[k] += eta[k]; + } } /* --- Store results -- -------------- */ atmos.backgrrecno[nspect] = mpi.backgrrecno; if (input.backgr_in_mem) { - storeBackground(nspect, 0, 0, chi_c, eta_c, sca_c); + storeBackground(nspect, 0, 0, chi_c, eta_c, sca_c); } else { - mpi.backgrrecno += writeBackground(nspect, 0, 0, - chi_c, eta_c, sca_c, NULL); + mpi.backgrrecno += writeBackground(nspect, 0, 0, + chi_c, eta_c, sca_c, NULL); } } } diff --git a/rh15d/iterate_p.c b/rh15d/iterate_p.c index 3c96c1a..705dd5b 100644 --- a/rh15d/iterate_p.c +++ b/rh15d/iterate_p.c @@ -141,7 +141,7 @@ void Iterate_p(int NmaxIter, double iterLimit) niter++; if (input.solve_ne == ITERATION) - Background(write_analyze_output=TRUE, equilibria_only=FALSE); + Background_p(write_analyze_output=TRUE, equilibria_only=FALSE); /* Update collisional radiative switching */ if (input.crsw > 0) diff --git a/rh15d/parallel.c b/rh15d/parallel.c index 7e54591..e2e36e3 100644 --- a/rh15d/parallel.c +++ b/rh15d/parallel.c @@ -112,6 +112,7 @@ void initParallelIO(bool_t run_ray, bool_t writej) { /* ------- begin -------------------------- closeParallelIO.c --- */ void closeParallelIO(bool_t run_ray, bool_t writej) { + int i; if (!run_ray) { close_hdf5_indata(); @@ -128,6 +129,12 @@ void closeParallelIO(bool_t run_ray, bool_t writej) { /* Free buffer variables */ freeBufVars(writej); + + /* Free RLK stuff */ + for (i = 0; i < atmos.Nrlk; i++) { + if (atmos.rlk_lines[i].zm != NULL) freeZeeman(atmos.rlk_lines[i].zm); + } + if (atmos.rlk_lines != NULL) free(atmos.rlk_lines); } /* ------- end -------------------------- closeParallelIO.c --- */ diff --git a/rh15d/writeRay.c b/rh15d/writeRay.c index 34cc622..0bd2fd6 100644 --- a/rh15d/writeRay.c +++ b/rh15d/writeRay.c @@ -656,14 +656,6 @@ void calculate_ray(void) { for (nact = 0; nact < atmos.Nactiveatom; nact++) { atom = atmos.activeatoms[nact]; - // TIAGO: commented below - /* Rewind atom files to point just before collisional data - if ((ierror = fseek(atom->fp_input, io.atom_file_pos[nact], SEEK_SET))) { - sprintf(messageStr, "Unable to rewind atom file for %s", atom->ID); - Error(ERROR_LEVEL_2, "rh15d_ray", messageStr); - } - */ - /* Free collisions matrix, going to be re-allocated in background */ if (atom->C != NULL) { freeMatrix((void **) atom->C); @@ -743,9 +735,7 @@ void calculate_ray(void) { if (input.prdh_limit_mem) prdh_limit_mem_save = TRUE; input.prdh_limit_mem = TRUE; } - Background_p(analyze_output=FALSE, equilibria_only=FALSE); - /* --- Solve radiative transfer for ray -- -------------- */ solveSpectrum(FALSE, FALSE); diff --git a/zeeman.c b/zeeman.c index 565c89b..0f953b5 100644 --- a/zeeman.c +++ b/zeeman.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Tue Jan 12 10:17:35 2010 -- + Last modified: Mon May 17 17:29:45 2021 -- -------------------------- ----------RH-- */ @@ -67,6 +67,7 @@ bool_t determinate(char *label, double g, int *n, double *S, int *L, *S = (multiplicity - 1) / 2.0; /* --- Orbital quantum number -- -------------- */ + *L = getOrbital(orbit[0]); /* --- Total angular momentum -- -------------- */ @@ -125,7 +126,7 @@ double effectiveLande(AtomicLine *line) g_l = Lande(S_l, L_l, J_l); g_u = Lande(S_u, L_u, J_u); - return 0.5*(g_u + g_l) + 0.25*(g_u - g_l) * + return 0.5*(g_u + g_l) + 0.25*(g_u - g_l) * (J_u*(J_u + 1.0) - J_l*(J_l + 1.0)); } else return 0.0; @@ -136,10 +137,7 @@ double effectiveLande(AtomicLine *line) double Lande(double S, int L, double J) { - if (J == 0.0) - return 0.0; - else - return 1.5 + (S*(S + 1.0) - L*(L + 1)) / (2.0*J*(J + 1.0)); + return 1.0 + zm_gamma(J, S, L); } /* ------- end ---------------------------- Lande.c ----------------- */ @@ -150,7 +148,7 @@ double ZeemanStrength(double Ju, double Mu, double Jl, double Ml) const char routineName[] = "ZeemanStrength"; int q, dJ; - double s; + double s = 0.0; /* --- Return the strength of Zeeman component (Ju, Mu) -> (Jl, Ml), where J and M are the total angular momentum and magnetic @@ -194,15 +192,17 @@ double ZeemanStrength(double Ju, double Mu, double Jl, double Ml) default: sprintf(messageStr, "Invalid dJ: %d", dJ); Error(ERROR_LEVEL_2, routineName, messageStr); - } + } return s; } /* ------- end ---------------------------- ZeemanStrength.c -------- */ /* ------- begin -------------------------- Zeeman.c ---------------- */ -ZeemanMultiplet* Zeeman(AtomicLine *line) +void Zeeman(AtomicLine *line) { + const char routineName[] = "Zeeman"; + register int n; bool_t result = TRUE; @@ -210,14 +210,14 @@ ZeemanMultiplet* Zeeman(AtomicLine *line) double Sl, Su, Jl, Ju, Mu, Ml, norm[3], gLu, gLl; Atom *atom = line->atom; ZeemanMultiplet *zm; - + /* --- Return a pointer to a ZeemanMultiplet structure with all the components of a Zeeman split line. The strengths in the line are normalized to unity for each of the three possible values of q = [-1, 0, 1]. Convention: - + -- q = +1 corresponds to a redshifted \sigma profile (zm->shift > 0). This redshifted profile has right-handed circular polarization when the @@ -228,7 +228,9 @@ ZeemanMultiplet* Zeeman(AtomicLine *line) -- -------------- */ - zm = (ZeemanMultiplet *) malloc(sizeof(ZeemanMultiplet)); + line->zm = (ZeemanMultiplet *) malloc(sizeof(ZeemanMultiplet)); + initZeeman(line->zm); + zm = line->zm; if (line->g_Lande_eff != 0.0) { @@ -282,21 +284,33 @@ ZeemanMultiplet* Zeeman(AtomicLine *line) if (fabs(Mu - Ml) <= 1.0) { zm->q[n] = (int) (Ml - Mu); zm->shift[n] = gLl*Ml - gLu*Mu; - zm->strength[n] = ZeemanStrength(Ju, Mu, Jl, Ml); - + zm->strength[n] = ZeemanStrength(Ju, Mu, Jl, Ml); + norm[zm->q[n]+1] += zm->strength[n]; - n++; + n++; } } } for (n = 0; n < zm->Ncomponent; n++) zm->strength[n] /= norm[zm->q[n]+1]; } - - return zm; + zm->g_eff = effectiveLande(line); } /* ------- end ---------------------------- Zeeman.c ---------------- */ +/* ------- begin -------------------------- zm_gamma.c -------------- */ + +double zm_gamma(double J, double S, double L) +{ + if (J == 0.0) + return 0.0; + else + return (J * (J + 1.0) + S * (S + 1.0) - L * (L + 1.0)) / + (2.0*J * (J + 1.0)); +} + +/* ------- end ---------------------------- zm_gamma.c -------------- */ + /* ------- begin -------------------------- adjustStokesMode.c ------ */ void adjustStokesMode() @@ -308,12 +322,12 @@ void adjustStokesMode() AtomicLine *line; /* --- Reset Stokes mode so that the full Stokes equations are - solved in case of input.StokesMode == FIELD_FREE, or - input.StokesMode == POLARIZATION_FREE. -- -------------- */ + solved in case of input.StokesMode == FIELD_FREE -- + -------------- */ - if (atmos.Nactiveatom == 0 || atmos.Stokes == FALSE || - input.StokesMode == NO_STOKES || - input.StokesMode == FULL_STOKES) return; + if ((atmos.Nactiveatom == 0 || atmos.Stokes == FALSE || + input.StokesMode == NO_STOKES || + input.StokesMode == FULL_STOKES) && !input.backgr_pol) return; else input.StokesMode = FULL_STOKES; @@ -331,7 +345,7 @@ void adjustStokesMode() /* --- First free up the space used in field-free calculation -- -------------- */ - + if (!input.limit_memory) freeMatrix((void **) line->phi); free(line->wphi); @@ -344,6 +358,19 @@ void adjustStokesMode() } /* ------- end ---------------------------- adjustStokesMode.c ------ */ +/* ------- begin -------------------------- initZeeman.c ------------ */ + +void initZeeman(ZeemanMultiplet *zm) +{ + zm->Ncomponent = 0; + zm->q = NULL; + + zm->shift = NULL; + zm->strength = NULL; + zm->g_eff = 0.0; +} +/* ------- end ---------------------------- initZeeman.c ------------ */ + /* ------- begin -------------------------- freeZeeman.c ------------ */ void freeZeeman(ZeemanMultiplet *zm) @@ -362,31 +389,30 @@ int getOrbital(char orbit) { const char routineName[] = "getOrbital"; - int L; + int L = 0; switch (orbit) { - case 'S': L = 0; break; - case 'P': L = 1; break; - case 'D': L = 2; break; - case 'F': L = 3; break; - case 'G': L = 4; break; - case 'H': L = 5; break; - case 'I': L = 6; break; - case 'K': L = 7; break; - case 'L': L = 8; break; - case 'M': L = 9; break; - case 'N': L = 10; break; - case 'O': L = 11; break; - case 'Q': L = 12; break; - case 'R': L = 13; break; - case 'T': L = 14; break; - case 'U': L = 15; break; - case 'V': L = 16; break; - case 'W': L = 17; break; - case 'X': L = 18; break; - case 'Y': L = 19; break; - case 'Z': L = 20; break; - default: + case 'S': L = 0; break; + case 'P': L = 1; break; + case 'D': L = 2; break; + case 'F': L = 3; break; + case 'G': L = 4; break; + case 'H': L = 5; break; + case 'I': L = 6; break; + case 'J': L = 7; break; + case 'K': L = 8; break; + case 'L': L = 9; break; + case 'M': L = 10; break; + case 'N': L = 11; break; + case 'O': L = 12; break; + case 'Q': L = 13; break; + case 'R': L = 14; break; + case 'T': L = 15; break; + case 'U': L = 16; break; + case 'V': L = 17; break; + case 'W': L = 18; break; + case 'X': L = 19; break; + default: sprintf(messageStr, "Invalid orbital: %c", orbit); Error(ERROR_LEVEL_2, routineName, messageStr); }