From 0226ac06ae0ce9cd6bdd263419ca4472ec0587b0 Mon Sep 17 00:00:00 2001 From: Tiago Pereira Date: Tue, 31 Jan 2012 19:45:55 +0000 Subject: [PATCH] background.c, chemequil.c, fpehandler.c, giigen.c, hydrogen.c, initial_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. --- atmos.h | 4 +- atom.h | 5 +- background.c | 69 +++- chemequil.c | 20 +- constant.h | 7 +- fpehandler.c | 194 ++-------- giigen.c | 4 +- hydrogen.c | 4 +- initial_xdr.c | 5 +- inputs.h | 10 +- iterate.c | 7 +- kurucz.c | 48 ++- ltepops.c | 86 ++--- molzeeman.c | 19 +- planck.c | 5 +- readatom.c | 10 +- readinput.c | 8 +- readj.c | 6 +- readmolecule.c | 87 ++++- readvalue.c | 27 +- scatter.c | 971 +++++++++++++++++++++++++------------------------ 21 files changed, 820 insertions(+), 776 deletions(-) diff --git a/atmos.h b/atmos.h index 75a6ea6..f65d483 100644 --- a/atmos.h +++ b/atmos.h @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Wed Mar 3 10:56:48 2010 -- + Last modified: Fri Jul 8 15:36:04 2011 -- -------------------------- ----------RH-- */ @@ -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); diff --git a/atom.h b/atom.h index d1b1ec2..2b5175e 100644 --- a/atom.h +++ b/atom.h @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Tue Aug 4 16:26:49 2009 -- + Last modified: Mon Apr 18 06:31:57 2011 -- -------------------------- ----------RH-- */ @@ -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}; @@ -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; diff --git a/background.c b/background.c index cca58b3..d7c6a4c 100644 --- a/background.c +++ b/background.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Fri Mar 5 09:16:05 2010 -- + Last modified: Thu Sep 29 21:36:56 2011 -- -------------------------- ----------RH-- */ @@ -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 @@ -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; @@ -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) @@ -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)); @@ -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 -- -------------- */ @@ -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 -- ------- */ @@ -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) -- */ @@ -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)) { @@ -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 -- -------------- */ @@ -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) { @@ -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) { @@ -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 -- -------------- */ @@ -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 -- -------------- */ @@ -592,7 +615,7 @@ void Background(bool_t analyzeoutput, bool_t equilibria_only) } } - if (analyzeoutput) { + if (write_analyze_output) { /* --- Write background record structure -- ------------ */ writeBRS(); @@ -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")) { @@ -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); diff --git a/chemequil.c b/chemequil.c index 678e564..d21ea6a 100644 --- a/chemequil.c +++ b/chemequil.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Tue Mar 31 09:15:25 2009 -- + Last modified: Mon Apr 18 07:08:41 2011 -- -------------------------- ----------RH-- */ @@ -23,9 +23,6 @@ #define COMMENT_CHAR "#" -/* --- Ionization energy Hmin in [J] -- -------------- */ - -#define E_ION_HMIN 0.754 * EV /* --- Acceleration parameters -- -------------- */ @@ -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 -- -------------- */ @@ -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: @@ -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); diff --git a/constant.h b/constant.h index e36ecab..1d42c8f 100644 --- a/constant.h +++ b/constant.h @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Wed Sep 14 22:13:01 2005 -- + Last modified: Wed Nov 17 16:27:27 2010 -- -------------------------- ----------RH-- */ @@ -60,6 +60,7 @@ #endif #define SQRTPI 1.77245385090551 + /* --- 1/(2sqrt(2)), needed for anisotropy of radiation -- ---------- */ #define TWOSQRTTWO 0.35355339059327 @@ -67,6 +68,10 @@ #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 -------------- */ diff --git a/fpehandler.c b/fpehandler.c index 2618a50..e6fd6b1 100644 --- a/fpehandler.c +++ b/fpehandler.c @@ -2,18 +2,12 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Tue Oct 15 16:26:36 2002 -- + Last modified: Mon Nov 22 15:09:15 2010 -- -------------------------- ----------RH-- */ /* --- Trap floating point exceptions on various machines: - SGI mips, IRIX 5.2 (compile with: -Dmips -DIRIX5) - SUN sparc, SunOS 4.x (compile with: -Dsparc -DSunOS4) - SUN sparc, SunOS 5.x (compile with: -Dsparc -DSunOS5) - DEC alpha, OSF1 3.x (compile with: -Dalpha -DOSF1V3) - DEC alpha, OSF1 3.x (compile with: -Dalpha -DOSF1V4) - -- -------------- */ #include @@ -22,112 +16,36 @@ #include "rh.h" #include "error.h" -extern char messageStr[]; - -#if !defined(SETNOTRAPS) - -#if defined(mips) +extern char messageStr[]; -/* --- Traps floating point exceptions for mips/IRIX5. - - CPU: mips - OS: IRIX5 - - see: ``man handle_sigfpes'' - - Requires linking with -lfpe - -- -------------- */ - -#include -#include +#if defined(SETNOTRAPS) +/* --- If SETNOTRAPS has been defined (see Makefile) -- ------------- */ void SetFPEtraps(void) { - sigfpe_[_UNDERFL].repls = _ZERO; - - sigfpe_[_OVERFL].abort = sigfpe_[_OVERFL].trace = 1; - sigfpe_[_DIVZERO].abort = sigfpe_[_DIVZERO].trace = 1; - sigfpe_[_INVALID].abort = sigfpe_[_INVALID].trace = 1; + /* --- Explicitly do not set traps -- -------------- */ - handle_sigfpes(_ON, _EN_OVERFL | _EN_DIVZERO | _EN_INVALID, - 0, _ABORT_ON_ERROR, 0); Error(MESSAGE, "SetFPEtraps", - "\n-Setting FPE traps for mips (IRIX 5.x)\n"); -} - -/* ------- end ------------------------------------------------------ */ - -#elif defined(sparc) - -#include - -#if defined(SunOS4) - -/* --- Traps floating point exceptions for sparc/SunOS4. - - CPU: sparc - OS: SunOS4 - - see: ``man ieee_handler'' - - -- -------------- */ - -void Trapped_FPE_Exception(int sig, int code, struct sigcontext *scp, - char *address); -int ieee_handler(const char *action, const char *exception, - sigfpe_handler_type hdl); - -void SetFPEtraps(void) -{ - const char routineName[] = "SetFPEtraps"; - - if (ieee_handler("set", "common", - (sigfpe_handler_type) Trapped_FPE_Exception) != 0) - Error(MESSAGE, routineName, "IEEE trapping not supported here\n"); - else - Error(MESSAGE, routineName, - "\n-Setting FPE traps for sparc (SunOS 4.x)\n"); + "\nFPE traps have not been set explicitly\n"); } +/* --- end SETNOTRAPS -- -------------- */ -void Trapped_FPE_Exception(int sig, int code, struct sigcontext *scp, - char *address) -{ - const char routineName[] = "Trapped_FPE_Exception"; - char *type = "unknown"; - - switch(code) { - - case FPE_INTDIV_TRAP: type = "integer divide by zero "; break; - case FPE_INTOVF_TRAP: type = "integer overflow "; break; - - case FPE_FLTDIV_TRAP: type = "floating point division by zero "; break; - case FPE_FLTUND_TRAP: type = "floating point underflow "; break; - case FPE_FLTOVF_TRAP: type = "floating point overflow "; break; - case FPE_FLTINEX_TRAP: type = "floating point inexact "; break; - } - - sprintf(messageStr, " ---- trapped IEEE FPE: %s ----\n" - " ---- signal: %d, code: 0x%x, addr: 0x%x, aborting ----\n", - type, sig, code, address); - Error(MESSAGE, routineName, messageStr); - abort(); -} -/* ------- end ------------------------------------------------------ */ +#else +/* --- SunOS -- -------------- */ -#else +#if defined(SunOS) -/* --- Traps floating point exceptions for sparc/SunOS5. +/* --- Traps floating point exceptions for SunOS5. - CPU: sparc OS: SunOS5 see: ``man ieee_handler'' Requires linking with -lsunmath -- -------------- */ - + #include #include #include @@ -171,98 +89,42 @@ void Trapped_FPE_Exception(int sig, siginfo_t *sip, ucontext_t *uap) Error(MESSAGE, routineName, messageStr); abort(); } -#endif - -/* ------- end ------------------------------------------------------ */ - - -#elif defined(alpha) -/* --- Traps floating point exceptions for DEC alpha/OSF1V[34] +/* --- end SunOS -- -------------- */ - CPU: alpha - OS: OSF1V[34] +/* --- Linux -- -------------- */ - see: ``man ieee'' +#elif defined(Linux) - -- -------------- */ - -#include -#include +#define _GNU_SOURCE 1 +#include -void Trapped_FPE_Exception(int dummy); -void SetFPEtraps() +void SetFPEtraps(void) { - const char routineName[] = "SetFPEtraps"; - - unsigned long fp_control = 0; - struct sigaction action, o_action; - - fp_control = (IEEE_TRAP_ENABLE_INV | - IEEE_TRAP_ENABLE_DZE | - IEEE_TRAP_ENABLE_OVF | - IEEE_TRAP_ENABLE_UNF); - ieee_set_fp_control(fp_control); - action.sa_handler = Trapped_FPE_Exception; + /* --- Enable some exceptions. + At startup all exceptions are masked. -- -------------- */ - if (sigaction(SIGFPE, &action, &o_action) != 0) - Error(MESSAGE, routineName, "\nIEEE trapping not supported here\n"); - else - Error(MESSAGE, routineName, - "\nSetting FPE traps for DEC alpha (OSF1 3.x)\n"); + feenableexcept(FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW); } -void Trapped_FPE_Exception(int dummy) -{ - const char routineName[] = "Trapped_FPE_Exception"; - - unsigned long fp_control = ieee_get_fp_control(); +/* --- end Linux -- -------------- */ - sprintf(messageStr, " ---- fp_control = 0x%lx\n", fp_control); - Error(MESSAGE, messageStr); - - if (fp_control & IEEE_STATUS_INV) - Error(MESSAGE, routineName, - " ---- trapped IEEE FPE: Invalid operation\n"); - if (fp_control & IEEE_STATUS_DZE) - Error(MESSAGE, routineName, - " ---- trapped IEEE FPE: floating point division by zero\n"); - if (fp_control & IEEE_STATUS_INE) - Error(MESSAGE, routineName, - " ---- trapped IEEE FPE: floating point inexact\n"); - if (fp_control & IEEE_STATUS_UNF) - Error(MESSAGE, routineName, - " ---- trapped IEEE FPE: floating point underflow\n"); - if (fp_control & IEEE_STATUS_OVF) - Error(MESSAGE, routineName, - " ---- trapped IEEE FPE: floating point overflow\n"); - - abort(); -} -/* ------- end ------------------------------------------------------ */ - -#else /* if defined(mips, sparc, alpha) */ +#else -/* ------- unknown -------------------------------------------------- */ +/* --- unknown -- -------------- */ void SetFPEtraps(void) { +/* Error(MESSAGE, "SetFPEtraps", "\nUnsupported CPU and/or OS: cannot set FPE traps explicitly\n"); +*/ } #endif -#else - -void SetFPEtraps(void) -{ - /* --- Explicitly do not set traps -- -------------- */ +/* --- end else -- -------------- */ - Error(MESSAGE, "SetFPEtraps", - "\nFPE traps have not been set explicitly\n"); -} - -#endif /* #if !defined(SETNOTRAPS) */ +#endif /* ------- end ------------------------------------------------------ */ diff --git a/giigen.c b/giigen.c index 0756463..3e7a9c4 100644 --- a/giigen.c +++ b/giigen.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Fri Jun 29 10:43:25 2001 -- + Last modified: Thu Jul 1 05:11:13 2010 -- -------------------------- ----------RH-- */ @@ -115,7 +115,7 @@ double GII(double adamp, double waveratio, double q_emit, double q_abs) if (q_emit >= PRD_QCORE) { phicore = exp(-SQ(q_emit)); - phiwing = adamp / (SQRTPI*(SQ(adamp) + SQ(q_emit))); + phiwing = adamp / (SQRTPI * (SQ(adamp) + SQ(q_emit))); pcore = phicore / (phicore + phiwing); } } diff --git a/hydrogen.c b/hydrogen.c index 4ed93b7..0ee4534 100644 --- a/hydrogen.c +++ b/hydrogen.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Fri Sep 25 14:07:18 2009 -- + Last modified: Tue Nov 16 12:42:01 2010 -- -------------------------- ----------RH-- */ @@ -91,6 +91,7 @@ void distribute_nH() -- -------------- */ if (atmos.H_LTE) { + atmos.H->NLTEpops = FALSE; Error(MESSAGE, routineName, "\nUsing LTE hydrogen populations for background opacities\n\n"); @@ -106,6 +107,7 @@ void distribute_nH() } } } else { + atmos.H->NLTEpops = TRUE; if (!atmos.H->active) atmos.H->n = matrix_double(atmos.H->Nlevel, atmos.Nspace); diff --git a/initial_xdr.c b/initial_xdr.c index c09d710..66a0cd2 100644 --- a/initial_xdr.c +++ b/initial_xdr.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Tue Jul 21 12:31:49 2009 -- + Last modified: Mon Jan 16 17:29:31 2012 -- -------------------------- ----------RH-- */ @@ -260,7 +260,8 @@ void initSolution(Atom *atom, Molecule *molecule) index = 0; for (nspect = 0; nspect < spectrum.Nspect; nspect++) { if (containsPRDline(&spectrum.as[nspect])) { - spectrum.PRDindex[nspect] = index++; + spectrum.PRDindex[nspect] = index; + index++; } } } diff --git a/inputs.h b/inputs.h index 5eb8b09..77c0652 100644 --- a/inputs.h +++ b/inputs.h @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Fri Apr 17 13:40:21 2009 -- + Last modified: Fri Jul 8 14:57:57 2011 -- -------------------------- ----------RH-- */ @@ -23,6 +23,7 @@ enum keywordtype {KEYWORD_REQUIRED, KEYWORD_DEFAULT, KEYWORD_OPTIONAL}; enum order_3D {LINEAR_3D, BICUBIC_3D}; +enum ne_solution {NONE, ONCE, ITERATION}; typedef struct { @@ -70,11 +71,13 @@ typedef struct { dampingFile[MAX_VALUE_LENGTH], coolingFile[MAX_VALUE_LENGTH], Itop[MAX_VALUE_LENGTH]; - bool_t magneto_optical, XRD, Eddington, solve_ne, - backgr_pol, limit_memory, allow_passive_bb, NonICE; + bool_t magneto_optical, XRD, Eddington, + backgr_pol, limit_memory, allow_passive_bb, NonICE, + rlkscatter; enum solution startJ; enum StokesMode StokesMode; enum order_3D interpolate_3D; + enum ne_solution solve_ne; enum PRDangle PRD_angle_dep; int isum, Ngdelay, Ngorder, Ngperiod, NmaxIter, PRD_NmaxIter, PRD_Ngdelay, PRD_Ngorder, PRD_Ngperiod, @@ -111,6 +114,7 @@ void setboolValue(char *value, void *pointer); void setintValue(char *value, void *pointer); void setdoubleValue(char *value, void *pointer); void setstartValue(char *value, void *pointer); +void setnesolution(char *value, void *pointer); void setStokesMode(char *value, void *pointer); void setPRDangle(char *value, void *pointer); void setThreadValue(char *value, void *pointer); diff --git a/iterate.c b/iterate.c index b9e74b5..8293710 100644 --- a/iterate.c +++ b/iterate.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Wed Apr 1 14:08:31 2009 -- + Last modified: Tue Nov 16 15:31:48 2010 -- -------------------------- ----------RH-- */ @@ -51,7 +51,7 @@ void Iterate(int NmaxIter, double iterLimit) register int niter, nact; double cswitch; - bool_t eval_operator; + bool_t eval_operator, write_analyze_output, equilibria_only; double dpopsmax, PRDiterlimit; Atom *atom; Molecule *molecule; @@ -128,6 +128,9 @@ void Iterate(int NmaxIter, double iterLimit) if (dpopsmax < iterLimit && cswitch <= 1.0 && input.prdswitch == 1.0 ) break; niter++; + if (input.solve_ne == ITERATION) + Background(write_analyze_output=TRUE, equilibria_only=FALSE); + /* Update collisional multiplier factor */ if (input.crsw > 0) cswitch = MAX(1.0, cswitch * pow(0.1, 1./input.crsw)); diff --git a/kurucz.c b/kurucz.c index dfca1f0..1feafe3 100644 --- a/kurucz.c +++ b/kurucz.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Wed Sep 30 16:04:28 2009 -- + Last modified: Tue Jul 12 10:57:37 2011 -- -------------------------- ----------RH-- */ @@ -412,7 +412,7 @@ void rlk_locate(int N, RLK_Line *lines, double lambda, int *low) /* ------- begin -------------------------- rlk_opacity.c ----------- */ 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) { register int k, n, kr; @@ -422,7 +422,8 @@ flags rlk_opacity(double lambda, int nspect, int mu, bool_t to_obs, 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, - phi, phi_Q, phi_U, phi_V, psi_Q, psi_U, psi_V; + phi, phi_Q, phi_U, phi_V, psi_Q, psi_U, psi_V, + epsilon, C, C2_atom, C2_ion, C3, dE, x; Atom *metal; AtomicLine *line; Element *element; @@ -447,6 +448,13 @@ flags rlk_opacity(double lambda, int nspect, int mu, bool_t to_obs, fourPI = 4.0 * PI; hc_4PI = hc / fourPI; + if (input.rlkscatter) { + C = 2 * PI * (Q_ELECTRON/EPSILON_0) * + (Q_ELECTRON/M_ELECTRON) / CLIGHT; + C2_atom = 2.15E-6; + C2_ion = 3.96E-6; + } + pf = (double *) malloc(atmos.Nspace * sizeof(double)); /* --- locate wavelength lambda in table of lines -- -------------- */ @@ -490,6 +498,9 @@ flags rlk_opacity(double lambda, int nspect, int mu, bool_t to_obs, chi[k] = 0.0; eta[k] = 0.0; } + if (input.rlkscatter) { + for (k = 0; k < atmos.Nspace; k++) scatt[k] = 0.0; + } } /* --- Add opacities from lines at this wavelength -- ------------- */ @@ -530,6 +541,17 @@ flags rlk_opacity(double lambda, int nspect, int mu, bool_t to_obs, 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; + } /* --- Set flag that line is present at this wavelength -- -- */ backgrflags.hasline = TRUE; @@ -558,6 +580,16 @@ flags rlk_opacity(double lambda, int nspect, int mu, bool_t to_obs, 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))); + + scatt[k] += (1.0 - epsilon) * chi_l * phi; + chi_l *= epsilon; + eta_l *= epsilon; + } + chi[k] += chi_l * phi; eta[k] += eta_l * phi; @@ -578,10 +610,6 @@ flags rlk_opacity(double lambda, int nspect, int mu, bool_t to_obs, } } } - /* - printf(" Line %f contributes at wavelength index %d\n", - rlk->lambda0, nspect); - */ } } } @@ -600,9 +628,9 @@ double RLKProfile(RLK_Line *rlk, int k, int mu, bool_t to_obs, { register int nz; - double v, phi_sm, phi_sp, phi_pi, psi_sm, psi_sp, psi_pi, adamp, - vB, H, F, sv, phi_sigma, phi_delta, sign, sin2_gamma, phi, - psi_sigma, psi_delta, vbroad, vtherm, GvdW, *np; + double v, phi_sm, phi_sp, phi_pi, psi_sm, psi_sp, psi_pi, adamp, + vB, H, F, sv, phi_sigma, phi_delta, sign, sin2_gamma, phi, + psi_sigma, psi_delta, vbroad, vtherm, GvdW, *np; Element *element; /* --- Returns the normalized profile for a Kurucz line diff --git a/ltepops.c b/ltepops.c index dd4659e..9531138 100644 --- a/ltepops.c +++ b/ltepops.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Thu Oct 18 14:47:37 2001 -- + Last modified: Wed Nov 17 10:07:34 2010 -- -------------------------- ----------RH-- */ @@ -37,7 +37,7 @@ void LTEpops(Atom *atom, bool_t Debeye) char labelStr[MAX_LINE_SIZE]; int Z, dZ, *nDebeye; long Nspace = atmos.Nspace; - double *cNe_T, *dE_kT, *dEion, dE, gi0, c1, *sum, c2; + double cNe_T, dE_kT, dEion, dE, gi0, c1, sum, c2; /* --- Computes LTE populations of a given atom. Takes account of Debeye shielding and lowering of the ionization @@ -69,55 +69,35 @@ void LTEpops(Atom *atom, bool_t Debeye) for (m = 1; m <= (atom->stage[i] - atom->stage[0]); m++, Z++) nDebeye[i] += Z; } - dEion = (double *) malloc(Nspace * sizeof(double)); - for (k = 0; k < Nspace; k++) - dEion[k] = c2 * sqrt(atmos.ne[k] / atmos.T[k]); - } else { - nDebeye = NULL; - dEion = NULL; } - - cNe_T = (double *) malloc(Nspace * sizeof(double)); - dE_kT = (double *) malloc(Nspace * sizeof(double)); - sum = (double *) malloc(Nspace * sizeof(double)); - /* --- Solve Saha-Boltzmann equilibrium equations -- ------------- */ for (k = 0; k < Nspace; k++) { - cNe_T[k] = 0.5*atmos.ne[k] * pow(c1/atmos.T[k], 1.5); - sum[k] = 1.0; - } + if (Debeye) dEion = c2 * sqrt(atmos.ne[k] / atmos.T[k]); + cNe_T = 0.5*atmos.ne[k] * pow(c1/atmos.T[k], 1.5); + sum = 1.0; - for (i = 1; i < atom->Nlevel; i++) { - dE = atom->E[i] - atom->E[0]; - gi0 = atom->g[i] / atom->g[0]; - dZ = atom->stage[i] - atom->stage[0]; - - if (Debeye) { - for (k = 0; k < Nspace; k++) - dE_kT[k] = (dE - nDebeye[i] * dEion[k]) / (KBOLTZMANN * atmos.T[k]); - } else { - for (k = 0; k < Nspace; k++) - dE_kT[k] = dE / (KBOLTZMANN * atmos.T[k]); - } - for (k = 0; k < Nspace; k++) { - atom->nstar[i][k] = gi0 * exp(-dE_kT[k]); - for (m = 1; m <= dZ; m++) atom->nstar[i][k] /= cNe_T[k]; - sum[k] += atom->nstar[i][k]; + for (i = 1; i < atom->Nlevel; i++) { + dE = atom->E[i] - atom->E[0]; + gi0 = atom->g[i] / atom->g[0]; + dZ = atom->stage[i] - atom->stage[0]; + + if (Debeye) + dE_kT = (dE - nDebeye[i] * dEion) / (KBOLTZMANN * atmos.T[k]); + else + dE_kT = dE / (KBOLTZMANN * atmos.T[k]); + + atom->nstar[i][k] = gi0 * exp(-dE_kT); + for (m = 1; m <= dZ; m++) atom->nstar[i][k] /= cNe_T; + sum += atom->nstar[i][k]; } - } - for (k = 0; k < Nspace; k++) { - atom->nstar[0][k] = atom->ntotal[k] / sum[k]; - } - for (i = 1; i < atom->Nlevel; i++) { - for (k = 0; k < Nspace; k++) - atom->nstar[i][k] *= atom->nstar[0][k]; - } + atom->nstar[0][k] = atom->ntotal[k] / sum; - if (Debeye) { - free(nDebeye); free(dEion); + for (i = 1; i < atom->Nlevel; i++) + atom->nstar[i][k] *= atom->nstar[0][k]; } - free(cNe_T); free(sum); free(dE_kT); + + if (Debeye) free(nDebeye); sprintf(labelStr, "LTEpops %2s", atom->ID); getCPU(3, TIME_POLL, labelStr); @@ -149,14 +129,16 @@ void LTEpops_elem(Element *element) } Linear(atmos.Npf, atmos.Tpf, element->pf[0], - atmos.Nspace, atmos.T, Uk, hunt=TRUE); + atmos.Nspace, atmos.T, Uk, hunt=TRUE); + for (i = 1; i < element->Nstage; i++) { Linear(atmos.Npf, atmos.Tpf, element->pf[i], - atmos.Nspace, atmos.T, Ukp1, hunt=TRUE); + atmos.Nspace, atmos.T, Ukp1, hunt=TRUE); for (k = 0; k < atmos.Nspace; k++) { element->n[i][k] = element->n[i-1][k] * CT_ne[k] * - exp(Ukp1[k] - Uk[k] - element->ionpot[i-1]/(KBOLTZMANN*atmos.T[k])); + exp(Ukp1[k] - Uk[k] - + element->ionpot[i-1]/(KBOLTZMANN*atmos.T[k])); sum[k] += element->n[i][k]; } SWAPPOINTER(Uk, Ukp1); @@ -190,7 +172,7 @@ void LTEmolecule(Molecule *molecule) register int k, v, J, kr; char labelStr[MAX_LINE_SIZE]; - double *kT, gJ, **E; + double kT, gJ, **E; MolecularLine *mrt; if (!molecule->active) { @@ -207,17 +189,15 @@ void LTEmolecule(Molecule *molecule) E[mrt->vj][(int) (mrt->gj - 1)/2] = mrt->Ej; } - kT = (double *) malloc(atmos.Nspace * sizeof(double)); - for (k = 0; k < atmos.Nspace; k++) { + for (k = 0; k < atmos.Nspace; k++) molecule->pf[k] = 0.0; - kT[k] = 1.0 / (KBOLTZMANN * atmos.T[k]); - } for (v = 0; v < molecule->Nv; v++) { for (J = 0; J < molecule->NJ; J++) { gJ = 2*J + 1; for (k = 0; k < atmos.Nspace; k++) - molecule->pfv[v][k] += gJ * exp(-E[v][J] * kT[k]); + molecule->pfv[v][k] += + gJ * exp(-E[v][J] / (KBOLTZMANN * atmos.T[k])); } /* --- Also store the total partition function here -- --------- */ @@ -225,7 +205,6 @@ void LTEmolecule(Molecule *molecule) molecule->pf[k] += molecule->pfv[v][k]; } - free(kT); freeMatrix((void **) E); sprintf(labelStr, "LTEpops %3s", molecule->ID); @@ -256,7 +235,6 @@ void SetLTEQuantities(void) atom. -- -------------- */ CollisionRate(atom, atom->fp_input); - fclose(atom->fp_input); /* --- Compute the fixed rates and store in Cij -- ------------ */ diff --git a/molzeeman.c b/molzeeman.c index c5ca6d3..9392b2e 100644 --- a/molzeeman.c +++ b/molzeeman.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Wed Apr 23 09:41:32 2003 -- + Last modified: Thu Jun 30 15:35:10 2011 -- -------------------------- ----------RH-- */ @@ -108,7 +108,7 @@ double MolLande_a(double Lambda, double Sigma, double Omega, double J) { /* --- Lande factor for level in Hund's case a - See: S.V. Berdyugina and S.K. Solanki 2002, A&A 385, 701-715, Eq. 10 + See: S.V. Berdyugina and S.K. Solanki 2002, A&A 385, 701-715, Eq. 2 -- -------------- */ return (Lambda + 2.0*Sigma) * Omega / (J * (J + 1.0)); @@ -121,13 +121,18 @@ double MolLande_b(double Lambda, double S, double N, double J) { /* --- Lande factor for level in Hund's case b - See: S.V. Berdyugina and S.K. Solanki 2002, A&A 385, 701-715, Eq. 2 + See: S.V. Berdyugina and S.K. Solanki 2002, A&A 385, 701-715, Eq. 10 -- -------------- */ - return 1.0 / (J*(J + 1.0)) * - (Lambda*Lambda / (2*N*(N + 1.0)) * - (J*(J + 1.0) + N*(N + 1.0) - S*(S + 1.0)) + - J*(J + 1.0) - N*(N + 1.0) + S*(S + 1.0)); + if (Lambda == 0) { + return 1.0 / (J*(J + 1.0)) * + (J*(J + 1.0) - N*(N + 1.0) + S*(S + 1.0)); + } else { + return 1.0 / (J*(J + 1.0)) * + (Lambda*Lambda / (2*N*(N + 1.0)) * + (J*(J + 1.0) + N*(N + 1.0) - S*(S + 1.0)) + + J*(J + 1.0) - N*(N + 1.0) + S*(S + 1.0)); + } } /* ------- end ---------------------------- MolLande_b.c ------------ */ diff --git a/planck.c b/planck.c index 3747c2c..8ea4670 100644 --- a/planck.c +++ b/planck.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Thu Feb 4 22:10:59 1999 -- + Last modified: Thu Nov 11 15:12:13 2010 -- -------------------------- ----------RH-- */ @@ -20,8 +20,11 @@ #include #include "rh.h" +#include "atom.h" +#include "spectrum.h" #include "constant.h" + #define MAX_EXPONENT 150.0 /* --- Function prototypes -- -------------- */ diff --git a/readatom.c b/readatom.c index 445d894..607e81b 100644 --- a/readatom.c +++ b/readatom.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Thu Jul 23 14:01:54 2009 -- + Last modified: Wed Nov 17 10:17:49 2010 -- -------------------------- ----------RH-- */ @@ -354,7 +354,6 @@ void readAtom(Atom *atom, char *atom_file, bool_t active) Nread = sscanf(inputLine, "%d %d %lf %d %s %lf", &j, &i, &continuum->alpha0, &continuum->Nlambda, nuDepStr, &lambdamin); - checkNread(Nread, Nrequired=6, routineName, checkPoint=5); continuum->j = MAX(i, j); continuum->i = MIN(i, j); @@ -540,6 +539,13 @@ void readAtom(Atom *atom, char *atom_file, bool_t active) atom->rhth = (rhthread *) malloc(input.Nthreads * sizeof(rhthread)); + /* --- Store the offset to allow pointing back to the start of the + collisional data in the atomic input file, and allocate + space for rate coefficients -- ------------- */ + + atom->offset_coll = ftell(atom->fp_input); + atom->C = matrix_double(SQ(Nlevel), atmos.Nspace); + } else { if (atom->popsFile) { diff --git a/readinput.c b/readinput.c index 855bc77..c40d184 100644 --- a/readinput.c +++ b/readinput.c @@ -2,7 +2,7 @@ Version: rh2.0, 1-D plane-parallel Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Fri Apr 17 13:29:42 2009 -- + Last modified: Mon Jul 11 09:12:48 2011 -- -------------------------- ----------RH-- */ @@ -108,10 +108,12 @@ void readInput() setboolValue}, {"KURUCZ_DATA", "none", FALSE, KEYWORD_OPTIONAL, &input.KuruczData, setcharValue}, + {"RLK_SCATTER", "FALSE", FALSE, KEYWORD_DEFAULT, &input.rlkscatter, + setboolValue}, {"KURUCZ_PF_DATA", "../../Atoms/pf_Kurucz.input", FALSE, KEYWORD_REQUIRED, &input.pfData, setcharValue}, - {"SOLVE_NE", "FALSE", FALSE, KEYWORD_DEFAULT, &input.solve_ne, - setboolValue}, + {"SOLVE_NE", "NONE", FALSE, KEYWORD_DEFAULT, &input.solve_ne, + setnesolution}, {"OPACITY_FUDGE", "none", FALSE, KEYWORD_OPTIONAL, &input.fudgeData, setcharValue}, {"METALLICITY", "0.0", FALSE, KEYWORD_DEFAULT, &input.metallicity, diff --git a/readj.c b/readj.c index 989f8ea..08d6ab9 100644 --- a/readj.c +++ b/readj.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Thu Jul 23 15:52:59 2009 -- + Last modified: Mon Jan 16 18:31:58 2012 -- -------------------------- ----------RH-- */ @@ -145,7 +145,7 @@ void readImu(int nspect, int mu, bool_t to_obs, double *I) index = spectrum.PRDindex[nspect]; offset = 2*(index*atmos.Nrays + mu) * recordsize; - if (to_obs) offset =+ recordsize; + if (to_obs) offset += recordsize; result &= (pread(spectrum.fd_Imu, I, recordsize, offset) == recordsize); @@ -174,7 +174,7 @@ void writeImu(int nspect, int mu, bool_t to_obs, double *I) index = spectrum.PRDindex[nspect]; offset = 2*(index*atmos.Nrays + mu) * recordsize; - if (to_obs) offset =+ recordsize; + if (to_obs) offset += recordsize; result &= (pwrite(spectrum.fd_Imu, I, recordsize, offset) == recordsize); diff --git a/readmolecule.c b/readmolecule.c index 486c09f..b7157cf 100644 --- a/readmolecule.c +++ b/readmolecule.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Wed Mar 3 14:45:32 2010 -- + Last modified: Wed Oct 12 20:00:35 2011 -- -------------------------- ----------RH-- */ @@ -30,9 +30,22 @@ #define COMMENT_CHAR "#" + + +/* --- Special case: C^13 fraction for CO -- -------------- */ + #define C13_FRACTION 0.011 +/* --- Special case: TI fractions for TiO -- -------------- */ + +#define TI_46 0.08 +#define TI_47 0.0743 +#define TI_48 0.738 +#define TI_49 0.055 +#define TI_50 0.054 + + /* --- Function prototypes -- -------------- */ int countTokens(char *line, char *separator); @@ -162,6 +175,8 @@ void readMolecule(Molecule *molecule, char *fileName, bool_t active) molecule->fit = KURUCZ_85; else if (strstr(fitStr, "SAUVAL_TATUM_84")) molecule->fit = SAUVAL_TATUM_84; + else if (strstr(fitStr, "IRWIN_81")) + molecule->fit = IRWIN_81; else if (strstr(fitStr, "TSUJI_73")) molecule->fit = TSUJI_73; else { @@ -483,6 +498,20 @@ int stringcompare(const void *s1, const void *s2) parity E 16 -- Isotope ID (16 refers to O16 in this case) + + -- KURUCZ_TIO + + Molecular line list for TiO molecules + + 708.9002 -1.213 68.0 9920.304 67.0 24022.776822a04p1 b08m1 46 +| | | | | | | | | | | | | | | | +0 5 10 20 30 40 50 60 70 + + Meaning of the numbers as above. + + Parity translation p --> E, and m --> F, for the orbital angular + momentum pointing parallel or anti-parallel to the internuclear axis. + -- -------------- */ void readMolecularLines(struct Molecule *molecule, char *line_data) @@ -592,7 +621,8 @@ void readMolecularLines(struct Molecule *molecule, char *line_data) mrt->polarizable = FALSE; } else if (strstr(format_string, "KURUCZ_CD18") || - strstr(format_string, "KURUCZ_NEW")) { + strstr(format_string, "KURUCZ_NEW") || + strstr(format_string, "KURUCZ_TIO")) { C = 2 * PI * (Q_ELECTRON/EPSILON_0) * (Q_ELECTRON/M_ELECTRON) / CLIGHT; @@ -606,9 +636,18 @@ void readMolecularLines(struct Molecule *molecule, char *line_data) mrt->isotope_frac = 1.0; getLine(fp_lines, COMMENT_CHAR, inputLine, exit_on_EOF=TRUE); - Nread = sscanf(inputLine, "%lf %lf %lf %lf %lf %lf", - &lambda_air, &log_gf, - &mrt->gi, &mrt->Ei, &mrt->gj, &mrt->Ej); + + if (strstr(format_string, "KURUCZ_TIO")) { + Nread = sscanf(inputLine, + "%10lf %7lf %5lf %10lf %5lf %10lf", + &lambda_air, &log_gf, + &mrt->gi, &mrt->Ei, &mrt->gj, &mrt->Ej); + } else { + Nread = sscanf(inputLine, + "%10lf %7lf %5lf %10lf %5lf %11lf", + &lambda_air, &log_gf, + &mrt->gi, &mrt->Ei, &mrt->gj, &mrt->Ej); + } checkNread(Nread, Nrequired=6, routineName, checkPoint=3); mrt->Ei = (HPLANCK * CLIGHT) / CM_TO_M * fabs(mrt->Ei); @@ -674,6 +713,44 @@ void readMolecularLines(struct Molecule *molecule, char *line_data) Nread = sscanf(inputLine+57, "%1d", &mrt->subi); Nread = sscanf(inputLine+65, "%1d", &mrt->subj); + + } else if (strstr(format_string, "KURUCZ_TIO")) { + + /* --- Electronic configuration -- -------------- */ + + strncpy(mrt->configi, inputLine+50, 2); + strncpy(mrt->configj, inputLine+59, 2); + mrt->configi[2] = '\0'; + mrt->configj[2] = '\0'; + + /* --- Parity designation -- -------------- */ + + mrt->parityi[0] = ((inputLine[53] == 'p') ? 'E' : 'F'); + mrt->parityj[0] = ((inputLine[61] == 'p') ? 'E' : 'F'); + mrt->parityi[1] = '\0'; + mrt->parityj[1] = '\0'; + + /* --- Vibrational quantum numbers -- ------------ */ + + Nread = sscanf(inputLine+51, "%2d", &mrt->vi); + Nread = sscanf(inputLine+59, "%2d", &mrt->vj); + + /* --- Subbranch (F1, F2, ...., E1, E2, ...etc) -- ---------- */ + + Nread = sscanf(inputLine+54, "%1d", &mrt->subi); + Nread = sscanf(inputLine+62, "%1d", &mrt->subj); + + /* --- Isotope fraction for Ti -- ------------ */ + + Nread = sscanf(inputLine+66, "%2d", &isotope); + switch (isotope) { + case 46: mrt->isotope_frac = TI_46; break; + case 47: mrt->isotope_frac = TI_47; break; + case 48: mrt->isotope_frac = TI_48; break; + case 49: mrt->isotope_frac = TI_49; break; + case 50: mrt->isotope_frac = TI_50; break; + default: mrt->isotope_frac = 1.0; + } } /* --- read additional data for Zeeman polarization -- -------- */ diff --git a/readvalue.c b/readvalue.c index 632bcb2..0f52a96 100644 --- a/readvalue.c +++ b/readvalue.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Fri Jul 24 13:23:13 2009 -- + Last modified: Mon Nov 22 11:20:47 2010 -- -------------------------- ----------RH-- */ @@ -381,6 +381,7 @@ void setInterpolate_3D(char *value, void *pointer) memcpy(pointer, &order, sizeof(enum order_3D)); } /* ------- end ---------------------------- setInterpolate_3D.c ----- */ + /* ------- begin -------------------------- setstartValue.c --------- */ void setstartValue(char *value, void *pointer) @@ -402,3 +403,27 @@ void setstartValue(char *value, void *pointer) memcpy(pointer, &startvalue, sizeof(enum_t)); } /* ------- end ---------------------------- setstartValue.c --------- */ + +/* ------- begin -------------------------- setnesolution.c --------- */ + +void setnesolution(char *value, void *pointer) +{ + const char routineName[] = "setnesolution"; + + enum ne_solution nesolution; + + if (!strcmp(value, "NONE")) + nesolution = NONE; + else if (!strcmp(value, "ONCE")) + nesolution = ONCE; + else if (!strcmp(value, "ITERATION")) { + nesolution = ITERATION; + } else { + sprintf(messageStr, + "Invalid value for keyword SOLVE_NE: %s", value); + Error(ERROR_LEVEL_2, routineName, messageStr); + } + + memcpy(pointer, &nesolution, sizeof(enum_t)); +} +/* ------- end ---------------------------- setnesolution.c --------- */ diff --git a/scatter.c b/scatter.c index c5e941e..e1fa6f7 100644 --- a/scatter.c +++ b/scatter.c @@ -2,7 +2,7 @@ Version: rh2.0 Author: Han Uitenbroek (huitenbroek@nso.edu) - Last modified: Fri Apr 4 22:32:56 2003 -- + Last modified: Mon Jan 16 17:45:32 2012 -- -------------------------- ----------RH-- */ @@ -299,28 +299,27 @@ void PRDScatter(AtomicLine *PRDline, enum Interpolation representation) } /* ------- end ---------------------------- PRDScatter.c ------------ */ +/* ------- begin -------------------------- PRDAngleScatter.c ------- */ -/* ------- begin -------------------------- PRDAngleApproxScatter.c - */ - -void PRDAngleApproxScatter(AtomicLine *PRDline, enum Interpolation representation) +void PRDAngleScatter(AtomicLine *PRDline, + enum Interpolation representation) { const char routineName[] = "scatterIntegral"; - register int la, k, lap, kr, ip, kxrd, mu, to_obs; + register int la, k, lap, kr, ip, mu, mup; char filename[MAX_LINE_SIZE]; - bool_t hunt, initialize; - int Np, Nread, Nwrite, ij, Nsubordinate, Nrays = atmos.Nrays; - long idx, lamu; - double q_emit, q0, qN, *q_abs = NULL, *qp = NULL, *wq = NULL, - *qpp = NULL, *gii = NULL, *adamp, cDop, gnorm, *J = NULL, - Jbar, scatInt, *J_k = NULL, *Pj, gamma, waveratio, vlos, sign, - *rho_intp, *rho_stor, *lambda_s, **rho_new; + bool_t hunt, initialize, to_obs, to_obs_p; + int Np, Nread, Nwrite, ij, lamu; + double *v_emit, v0, vN, *v_abs = NULL, *vp = NULL, *wv = NULL, + *rii = NULL, *adamp, *Jbar, cDop, *RIInorm, *I = NULL, + **Imup, *Ik, *Pj, *gamma, **v_los, *phi_emit, wmup, *sv; Atom *atom; - AtomicLine *line, **XRD, *XRDline; + AtomicLine *line; AtomicContinuum *continuum; - /* --- This the static case, with Jorrit Leenaart's approximation */ - + /* --- Calculate the angle-dependent scattering integral when + angle-dependent scattering is requested. -- -------------- */ + atom = PRDline->atom; if (!PRDline->PRD) { @@ -328,11 +327,12 @@ void PRDAngleApproxScatter(AtomicLine *PRDline, enum Interpolation representatio PRDline->j, PRDline->i, atom->ID); Error(ERROR_LEVEL_2, routineName, messageStr); } - getCPU(3, TIME_START, NULL); - /* --- Open temporary file for storage of redistribution weights when - called for the first time -- -------------- */ + cDop = (NM_TO_M * PRDline->lambda0) / (4.0 * PI); + + /* --- When called for the first time open temporary file for + storage of redistribution weights -- ------------ */ initialize = FALSE; if (PRDline->fp_GII == NULL) { @@ -355,310 +355,275 @@ void PRDAngleApproxScatter(AtomicLine *PRDline, enum Interpolation representatio } } if (!initialize) rewind(PRDline->fp_GII); - - /* --- Set XRD line array -- -------------- */ - - if (input.XRD) - Nsubordinate = 1 + PRDline->Nxrd; - else - Nsubordinate = 1; - - XRD = (AtomicLine **) malloc(Nsubordinate * sizeof(AtomicLine*)); - XRD[0] = PRDline; - if (input.XRD) { - for (kxrd = 0; kxrd < PRDline->Nxrd; kxrd++) - XRD[kxrd + 1] = PRDline->xrd[kxrd]; - } - /* --- Initialize the emission profile ratio rho -- -------------- */ - for (la = 0; la < PRDline->Nlambda; la++) { - for (mu = 0; mu < Nrays; mu++) { - for (to_obs = 0; to_obs <= 1; to_obs++) { - lamu = 2*(Nrays*la + mu) + to_obs; - for (k = 0; k < atmos.Nspace; k++) { - PRDline->rho_prd[lamu][k] = 1.0; - } - } - } - } - Pj = (double *) malloc (atmos.Nspace * sizeof(double)); - adamp = (double *) malloc (atmos.Nspace * sizeof(double)); - cDop = (NM_TO_M * PRDline->lambda0) / (4.0 * PI); + /* --- Temporary storage space -- -------------- */ - for (k = 0; k < atmos.Nspace; k++) { - adamp[k] = (PRDline->Grad + PRDline->Qelast[k]) * cDop / atom->vbroad[k]; + Imup = matrix_double(PRDline->Nlambda, atmos.Nspace); + Ik = (double *) malloc(PRDline->Nlambda * sizeof(double)); + v_emit = (double *) malloc(atmos.Nspace * sizeof(double)); + v_abs = (double *) malloc(PRDline->Nlambda * sizeof(double)); + adamp = (double *) malloc(atmos.Nspace * sizeof(double)); + v_los = matrix_double(atmos.Nrays, atmos.Nspace); + gamma = (double *) malloc(atmos.Nspace * sizeof(double)); + Pj = (double *) malloc(atmos.Nspace * sizeof(double)); + Jbar = (double *) malloc(atmos.Nspace * sizeof(double)); + RIInorm = (double *) malloc(atmos.Nspace * sizeof(double)); + sv = (double *) malloc(atmos.Nspace * sizeof(double)); - /* --- Evaluate the total rate Pj out of the line's upper level - */ + /* --- Evaluate first the total rate Pj out of the line's upper level + and then the coherency fraction gamma -- -------------- */ - Pj[k] = PRDline->Qelast[k]; - for (ip = 0; ip < atom->Nlevel; ip++) { - ij = ip * atom->Nlevel + PRDline->j; + for (k = 0; k < atmos.Nspace; k++) + Pj[k] = PRDline->Qelast[k]; + for (ip = 0; ip < atom->Nlevel; ip++) { + ij = ip * atom->Nlevel + PRDline->j; + for (k = 0; k < atmos.Nspace; k++) Pj[k] += atom->C[ij][k]; - } - for (kr = 0; kr < atom->Nline; kr++) { - line = &atom->line[kr]; - if (line->j == PRDline->j) Pj[k] += line->Rji[k]; - if (line->i == PRDline->j) Pj[k] += line->Rij[k]; - } - for (kr = 0; kr < atom->Ncont; kr++) { - continuum = &atom->continuum[kr]; - if (continuum->j == PRDline->j) Pj[k] += continuum->Rji[k]; - if (continuum->i == PRDline->j) Pj[k] += continuum->Rij[k]; - } } - /* --- Loop over subordinate lines -- -------------- */ - - for (kxrd = 0; kxrd < Nsubordinate; kxrd++) { - XRDline = XRD[kxrd]; - waveratio = XRDline->lambda0 / PRDline->lambda0; + for (kr = 0; kr < atom->Nline; kr++) { + line = &atom->line[kr]; + if (line->j == PRDline->j) + for (k = 0; k < atmos.Nspace; k++) Pj[k] += line->Rji[k]; + if (line->i == PRDline->j) + for (k = 0; k < atmos.Nspace; k++) Pj[k] += line->Rij[k]; + } + for (kr = 0; kr < atom->Ncont; kr++) { + continuum = &atom->continuum[kr]; + if (continuum->j == PRDline->j) + for (k = 0; k < atmos.Nspace; k++) Pj[k] += continuum->Rji[k]; + if (continuum->i == PRDline->j) + for (k = 0; k < atmos.Nspace; k++) Pj[k] += continuum->Rij[k]; + } + for (k = 0; k < atmos.Nspace; k++) { + gamma[k] = atom->n[PRDline->i][k] / atom->n[PRDline->j][k] * + PRDline->Bij / Pj[k]; + } + /* --- Store line-of-sight velocity in Doppler units to avoid + having to recompute it for every wavelength -- ------------- */ - J_k = realloc(J_k, XRDline->Nlambda * sizeof(double)); - q_abs = realloc(q_abs, XRDline->Nlambda * sizeof(double)); + for (mu = 0; mu < atmos.Nrays; mu++) { + for (k = 0; k < atmos.Nspace; k++) + v_los[mu][k] = vproject(k, mu) / atom->vbroad[k]; + } + /* --- Depth-dependent damping parameter -- -------------- */ - /* --- Loop over all spatial locations -- -------------- */ + for (k = 0; k < atmos.Nspace; k++) { + Jbar[k] = PRDline->Rij[k] / PRDline->Bij; + sv[k] = 1.0 / (SQRTPI * atom->vbroad[k]); + adamp[k] = + (PRDline->Grad + PRDline->Qelast[k]) * cDop / atom->vbroad[k]; + } + /* --- Outer loop over emission wavelength -- -------------- */ - for (k = 0; k < atmos.Nspace; k++) { - gamma = atom->n[XRDline->i][k] / atom->n[PRDline->j][k] * - XRDline->Bij / Pj[k]; - Jbar = XRDline->Rij[k] / XRDline->Bij; + for (la = 0; la < PRDline->Nlambda; la++) { - /* --- Get local mean intensity and wavelength in Doppler units */ + /* --- Loop over emission angle (down and up) -- -------------- */ - for (la = 0; la < XRDline->Nlambda; la++) { - // Tiago: this is the line where we need to replace spectrum.J by spectrum.Jgas + for (mu = 0; mu < atmos.Nrays; mu++) { + for (to_obs = 0; to_obs <= 1; to_obs++) { + lamu = 2*(atmos.Nrays*la + mu) + to_obs; - J_k[la] = spectrum.Jgas[XRDline->Nblue + la][k]; - //J_k[la] = spectrum.J[XRDline->Nblue + la][k]; - q_abs[la] = (XRDline->lambda[la] - XRDline->lambda0) * CLIGHT / - (XRDline->lambda0 * atom->vbroad[k]); - } - switch (representation) { - case LINEAR: - break; - case SPLINE: - splineCoef(XRDline->Nlambda, q_abs, J_k); - break; - case EXP_SPLINE: - exp_splineCoef(XRDline->Nlambda, q_abs, J_k, TENSION); - break; - } - /* --- Outer wavelength loop over emission wavelengths -- ----- */ + if (atmos.moving || + (PRDline->polarizable && input.StokesMode > FIELD_FREE)) + phi_emit = PRDline->phi[lamu]; + else + phi_emit = PRDline->phi[la]; - for (la = 0; la < PRDline->Nlambda; la++) { - q_emit = (PRDline->lambda[la] - PRDline->lambda0) * CLIGHT / - (PRDline->lambda0 * atom->vbroad[k]); + for (k = 0; k < atmos.Nspace; k++) { + PRDline->rho_prd[lamu][k] = 0.0; + RIInorm[k] = 0.0; + } + /* --- Depth-dependent emission wavelengths for this + direction (in Doppler units). -- -------------- */ - /* --- Establish integration limits over absorption wavelength, - using only regions where the redistribution function is - non-zero. (See also function GII.) -- -------------- */ - - if (fabs(q_emit) < PRD_QCORE) { - q0 = -PRD_QWING; - qN = PRD_QWING; + if (to_obs) { + for (k = 0; k < atmos.Nspace; k++) + v_emit[k] = (PRDline->lambda[la] - PRDline->lambda0) * CLIGHT / + (atom->vbroad[k] * PRDline->lambda0) + v_los[mu][k]; } else { - if (fabs(q_emit) < PRD_QWING) { - if (q_emit > 0.0) { - q0 = -PRD_QWING; - qN = waveratio * (q_emit + PRD_QSPREAD); - } else { - q0 = waveratio * (q_emit - PRD_QSPREAD); - qN = PRD_QWING; - } - } else { - q0 = waveratio * (q_emit - PRD_QSPREAD); - qN = waveratio * (q_emit + PRD_QSPREAD); - } + for (k = 0; k < atmos.Nspace; k++) + v_emit[k] = (PRDline->lambda[la] - PRDline->lambda0) * CLIGHT / + (atom->vbroad[k] * PRDline->lambda0) - v_los[mu][k]; } - Np = (int) ((qN - q0) / PRD_DQ) + 1; - qp = (double *) realloc(qp, Np * sizeof(double)); - for (lap = 1, qp[0] = q0; lap < Np; lap++) - qp[lap] = qp[lap - 1] + PRD_DQ; - /* --- Fold interpolation of J for symmetric lines. -- ------ */ + /* --- Loop over absorption directions -- -------------- */ - if (XRDline->symmetric) { - qpp = (double *) realloc(qpp, Np * sizeof(double)); - for (lap = 0; lap < Np; lap++) qpp[lap] = fabs(qp[lap]); - } + for (mup = 0; mup < atmos.Nrays; mup++) { + wmup = 0.5 * atmos.wmu[mup]; + for (to_obs_p = 0; to_obs_p <= 1; to_obs_p++) { - /* --- Interpolate mean intensity onto fine grid. Choose linear, - spline or exponential spline interpolation -- -------- */ + /* --- Read specific intensity in this direction for all + wavelengths -- -------------- */ - J = (double *) realloc(J, Np * sizeof(double)); - switch (representation) { - case LINEAR: - Linear(XRDline->Nlambda, q_abs, J_k, Np, - (XRDline->symmetric) ? qpp : qp, J, hunt=TRUE); - break; - case SPLINE: - splineEval(Np, (XRDline->symmetric) ? qpp : qp, J, hunt=TRUE); - break; - case EXP_SPLINE: - exp_splineEval(Np, (XRDline->symmetric) ? qpp : qp, J, hunt=TRUE); - break; - } - /* --- Compute the redistribution weights -- -------------- */ + for (lap = 0; lap < PRDline->Nlambda; lap++) + readImu(PRDline->Nblue + lap, mup, to_obs_p, Imup[lap]); - gii = (double *) realloc(gii, 2*Nrays*Np * sizeof(double)); - //gii = (double *) realloc(gii, Np * sizeof(double)); - - if (initialize) { + /* --- Loop over space -- -------------- */ - /* --- Integration weights (See: Press et al. Numerical Recipes, - p. 107, eq. 4.1.12) -- -------------- */ + for (k = 0; k < atmos.Nspace; k++) { + for (lap = 0; lap < PRDline->Nlambda; lap++) { + Ik[lap] = Imup[lap][k]; - wq = (double *) realloc(wq, Np * sizeof(double)); - wq[0] = 5.0/12.0 * PRD_DQ; - wq[1] = 13.0/12.0 * PRD_DQ; - for (lap = 2; lap < Np-2; lap++) wq[lap] = PRD_DQ; - wq[Np-1] = 5.0/12.0 * PRD_DQ; - wq[Np-2] = 13.0/12.0 * PRD_DQ; + /* --- Array of absorption wavelengths in this + direction -- -------------- */ - // Tiago: adding angle loops here - for (mu = 0; mu < Nrays; mu++) { - for (to_obs = 0; to_obs <= 1; to_obs++) { - - sign = (to_obs) ? 1.0 : -1.0; - - //vlos = sign*spectrum.v_los[mu][k] / atom->vbroad[k]; - vlos = 0.; - - for (lap = 0; lap < Np; lap++) { - - idx = mu*2*Np + to_obs*Np + lap; /* index of gii */ + if (to_obs_p) + v_abs[lap] = + (PRDline->lambda[lap] - PRDline->lambda0) * CLIGHT / + (PRDline->lambda0 * atom->vbroad[k]) + v_los[mup][k]; + else + v_abs[lap] = + (PRDline->lambda[lap] - PRDline->lambda0) * CLIGHT / + (PRDline->lambda0 * atom->vbroad[k]) - v_los[mup][k]; + } + /* --- Setup spline interpolation coefficients for the + integration over absorption intensity -- ------- */ - gii[idx] = GII(adamp[k], waveratio, q_emit-vlos, qp[lap]) * wq[lap]; + switch (representation) { + case SPLINE: + splineCoef(PRDline->Nlambda, v_abs, Ik); + break; + case EXP_SPLINE: + exp_splineCoef(PRDline->Nlambda, v_abs, Ik, TENSION); + break; + case LINEAR: break; } - } - } + /* --- Establish integration limits over absorption + wavelength, using only regions where the + redistribution function is non-zero. -- -------- */ + + if (fabs(v_emit[k]) < PRD_QCORE) { + v0 = -PRD_QWING; + vN = PRD_QWING; + } else { + if (fabs(v_emit[k]) < PRD_QWING) { + if (v_emit[k] > 0.0) { + v0 = -PRD_QWING; + vN = v_emit[k] + PRD_QSPREAD; + } else { + v0 = v_emit[k] - PRD_QSPREAD; + vN = PRD_QWING; + } + } else { + v0 = v_emit[k] - PRD_QSPREAD; + vN = v_emit[k] + PRD_QSPREAD; + } + } + Np = (int) ((vN - v0) / PRD_DQ) + 1; + vp = (double *) realloc(vp, Np * sizeof(double)); + for (lap = 1, vp[0] = v0; lap < Np; lap++) + vp[lap] = vp[lap - 1] + PRD_DQ; - if ((Nwrite = - //fwrite(gii, sizeof(double), Np, PRDline->fp_GII)) != Np) { - fwrite(gii, sizeof(double), 2*Nrays*Np, PRDline->fp_GII)) != 2*Nrays*Np) { - sprintf(messageStr, + /* --- Interpolate specific intensity onto fine grid. + Choose linear, spline or exponential spline + interpolation -- -------------- */ + + I = (double *) realloc(I, Np * sizeof(double)); + switch (representation) { + case LINEAR: + Linear(PRDline->Nlambda, v_abs, Ik, Np, vp, I, hunt=TRUE); + break; + case SPLINE: + splineEval(Np, vp, I, hunt=TRUE); + break; + case EXP_SPLINE: + exp_splineEval(Np, vp, I, hunt=TRUE); + break; + } + /* --- Compute the redistribution weights -- ---------- */ + + rii = (double *) realloc(rii, Np * sizeof(double)); + if (initialize) { + + /* --- Integration weights (See: Press et al. + Numerical Recipes, p. 107, eq. 4.1.12) -- ---- */ + + wv = (double *) realloc(wv, Np * sizeof(double)); + wv[0] = 5.0/12.0 * PRD_DQ; + wv[1] = 13.0/12.0 * PRD_DQ; + for (lap = 2; lap < Np-2; lap++) wv[lap] = PRD_DQ; + wv[Np-1] = 5.0/12.0 * PRD_DQ; + wv[Np-2] = 13.0/12.0 * PRD_DQ; + + /* --- RII(x, mu, x', mu')/phi(x, mu) -- ------------ */ + + for (lap = 0; lap < Np; lap++) { + rii[lap] = RII(v_emit[k], vp[lap], adamp[k], mu, mup) * + (sv[k] / phi_emit[k]) * wv[lap] * wmup; + } + if ((Nwrite = fwrite(rii, sizeof(double), Np, + PRDline->fp_GII)) != Np) { + sprintf(messageStr, "Unable to write proper number of redistribution weights\n" " Wrote %d instead of %d.\n Line %d -> %d, la = %d, k = %d", - Nwrite, Np, XRDline->j, XRDline->i, la, k); - Error(ERROR_LEVEL_2, routineName, messageStr); + Nwrite, Np, PRDline->j, PRDline->i, la, k); + Error(ERROR_LEVEL_2, routineName, messageStr); + } + } else { + if ((Nread = fread(rii, sizeof(double), Np, + PRDline->fp_GII)) != Np) { + sprintf(messageStr, + "Unable to read proper number of redistribution weights\n" + " Read %d instead of %d.\n Line %d -> %d, la = %d, k = %d", + Nread, Np, PRDline->j, PRDline->i, la, k); + Error(ERROR_LEVEL_2, routineName, messageStr); + } + } + /* --- Inner wavelength loop doing actual wavelength + integration over absorption wavelengths -- ----- */ + + for (lap = 0; lap < Np; lap++) { + RIInorm[k] += rii[lap]; + PRDline->rho_prd[lamu][k] += I[lap] * rii[lap]; + } + } } - } else { - if ((Nread = - //fread(gii, sizeof(double), Np, PRDline->fp_GII)) != Np) { - fread(gii, sizeof(double), 2*Nrays*Np, PRDline->fp_GII)) != 2*Nrays*Np) { - sprintf(messageStr, - "Unable to read proper number of redistribution weights\n" - " Read %d instead of %d.\n Line %d -> %d, la = %d, k = %d", - Nread, Np, XRDline->j, XRDline->i, la, k); - Error(ERROR_LEVEL_2, routineName, messageStr); - } } - /* --- Inner wavelength loop doing actual wavelength integration - over absorption wavelengths -- -------------- */ - for (mu = 0; mu < Nrays; mu++) { - for (to_obs = 0; to_obs <= 1; to_obs++) { - - lamu = 2*(Nrays*la + mu) + to_obs; - //lamu = la; - - gnorm = 0.0; - scatInt = 0.0; - - // Integral in absorption wavelength: - for (lap = 0; lap < Np; lap++) { - idx = mu*2*Np + to_obs*Np + lap; /* index of gii */ - //idx = lap; - gnorm += gii[idx]; - scatInt += J[lap] * gii[idx]; - } - - // Tiago: put condition to prevent negative scatInt? - PRDline->rho_prd[lamu][k] += gamma*(scatInt/gnorm - Jbar); /* Integral over XRD lines*/ - } /* end of to_obs loop */ - } /* end of mu loop */ - - // Tiago: must interpolate rho_prd back into the lab frame. use Linear() - } /* End of emission wavelength loop */ - } /* End of depth loop */ - } /* End of subordinate line loop */ - - /* --- Tiago: shift rho back to lab rest frame -- -------------- */ - - rho_stor = (double *) calloc(PRDline->Nlambda, sizeof(double)); - rho_intp = (double *) calloc(PRDline->Nlambda, sizeof(double)); - lambda_s = (double *) calloc(PRDline->Nlambda, sizeof(double)); - rho_new = matrix_double(2*atmos.Nrays * PRDline->Nlambda, atmos.Nspace); - - /* - for (k = 0; k < atmos.Nspace; k++) { - for (mu = 0; mu < Nrays; mu++ ) { - for (to_obs = 0; to_obs <= 1; to_obs++) { - - // Grab rho[lambda] and shifted wavelengths - for (la = 0; la < PRDline->Nlambda; la++) { - lamu = 2*(Nrays*la + mu) + to_obs; - sign = (to_obs) ? 1.0 : -1.0; - rho_stor[la] = PRDline->rho_prd[lamu][k]; - lambda_s[la] = PRDline->lambda[la]*(1.-spectrum.v_los[mu][k]*sign/CLIGHT); - } - - // Linearly interpolate, equivalent to wavelength shift - Linear(PRDline->Nlambda, lambda_s, rho_stor, PRDline->Nlambda, - PRDline->lambda, rho_intp, hunt=TRUE); - - // Save back into rho - for (la = 0; la < PRDline->Nlambda; la++) { - lamu = 2*(Nrays*la + mu) + to_obs; - rho_new[lamu][k] = rho_intp[la]; - PRDline->rho_prd[lamu][k] = rho_intp[la]; + for (k = 0; k < atmos.Nspace; k++) { + PRDline->rho_prd[lamu][k] = 1.0 + + gamma[k] * (PRDline->rho_prd[lamu][k]/RIInorm[k] - Jbar[k]); } - } // direction loop - } // mu loop - - } // depth loop - */ - - - free(rho_stor); free(rho_intp); free(lambda_s); freeMatrix((void **) rho_new); - - /* --- Clean temporary variable space -- -------------- */ - - for (kxrd = 0; kxrd < Nsubordinate; kxrd++) { - if (XRD[kxrd]->symmetric) { - free(qpp); - break; + } } } - free(J_k); free(q_abs); free(gii); - free(qp); free(wq); free(J); - free(XRD); free(Pj); free(adamp); + /* --- Clean temporary variable space -- -------------- */ + + freeMatrix((void **) Imup); + freeMatrix((void **) v_los); + + free(Ik); free(v_abs); free(v_emit); free(adamp); + free(gamma); free(Pj); + free(wv); free(I); free(rii); + free(vp); free(Jbar); free(RIInorm); free(sv); sprintf(messageStr, "Scatter Int %5.1f", PRDline->lambda0); getCPU(3, TIME_POLL, messageStr); } -/* ------- end ---------------------------- PRDAngleApproxScatter.c - */ +/* ------- end ---------------------------- PRDAngleScatter.c ------- */ -/* ------- begin -------------------------- PRDAngleScatter.c ------- */ -void PRDAngleScatter(AtomicLine *PRDline, - enum Interpolation representation) +/* ------- begin -------------------------- PRDAngleApproxScatter.c - */ + +void PRDAngleApproxScatter(AtomicLine *PRDline, enum Interpolation representation) { const char routineName[] = "scatterIntegral"; - register int la, k, lap, kr, ip, mu, mup; + register int la, k, lap, kr, ip, kxrd, mu, to_obs; char filename[MAX_LINE_SIZE]; - bool_t hunt, initialize, to_obs, to_obs_p; - int Np, Nread, Nwrite, ij, lamu; - double *v_emit, v0, vN, *v_abs = NULL, *vp = NULL, *wv = NULL, - *rii = NULL, *adamp, *Jbar, cDop, *RIInorm, *I = NULL, - **Imup, *Ik, *Pj, *gamma, **v_los, *phi_emit, wmup, *sv; + bool_t hunt, initialize; + int Np, Nread, Nwrite, ij, Nsubordinate, Nrays = atmos.Nrays; + long idx, lamu; + double q_emit, q0, qN, *q_abs = NULL, *qp = NULL, *wq = NULL, + *qpp = NULL, *gii = NULL, *adamp, cDop, gnorm, *J = NULL, + Jbar, scatInt, *J_k = NULL, *Pj, gamma, waveratio, vlos, sign, + *rho_intp, *rho_stor, *lambda_s, **rho_new; Atom *atom; - AtomicLine *line; + AtomicLine *line, **XRD, *XRDline; AtomicContinuum *continuum; - /* --- Calculate the angle-dependent scattering integral when - angle-dependent scattering is requested. -- -------------- */ - + /* --- This the static case, with Jorrit Leenaart's approximation */ + atom = PRDline->atom; if (!PRDline->PRD) { @@ -666,12 +631,11 @@ void PRDAngleScatter(AtomicLine *PRDline, PRDline->j, PRDline->i, atom->ID); Error(ERROR_LEVEL_2, routineName, messageStr); } + getCPU(3, TIME_START, NULL); - cDop = (NM_TO_M * PRDline->lambda0) / (4.0 * PI); - - /* --- When called for the first time open temporary file for - storage of redistribution weights -- ------------ */ + /* --- Open temporary file for storage of redistribution weights when + called for the first time -- -------------- */ initialize = FALSE; if (PRDline->fp_GII == NULL) { @@ -694,246 +658,285 @@ void PRDAngleScatter(AtomicLine *PRDline, } } if (!initialize) rewind(PRDline->fp_GII); + + /* --- Set XRD line array -- -------------- */ - /* --- Temporary storage space -- -------------- */ - - Imup = matrix_double(PRDline->Nlambda, atmos.Nspace); - Ik = (double *) malloc(PRDline->Nlambda * sizeof(double)); - v_emit = (double *) malloc(atmos.Nspace * sizeof(double)); - v_abs = (double *) malloc(PRDline->Nlambda * sizeof(double)); - adamp = (double *) malloc(atmos.Nspace * sizeof(double)); - v_los = matrix_double(atmos.Nrays, atmos.Nspace); - gamma = (double *) malloc(atmos.Nspace * sizeof(double)); - Jbar = (double *) malloc(atmos.Nspace * sizeof(double)); - RIInorm = (double *) malloc(atmos.Nspace * sizeof(double)); - sv = (double *) malloc(atmos.Nspace * sizeof(double)); - - /* --- Evaluate first the total rate Pj out of the line's upper level - and then the coherency fraction gamma -- -------------- */ + if (input.XRD) + Nsubordinate = 1 + PRDline->Nxrd; + else + Nsubordinate = 1; - Pj = gamma; - for (ip = 0; ip < atom->Nlevel; ip++) { - ij = ip * atom->Nlevel + PRDline->j; - for (k = 0; k < atmos.Nspace; k++) - Pj[k] = PRDline->Qelast[k] + atom->C[ij][k]; - } - for (kr = 0; kr < atom->Nline; kr++) { - line = &atom->line[kr]; - if (line->j == PRDline->j) - for (k = 0; k < atmos.Nspace; k++) Pj[k] += line->Rji[k]; - if (line->i == PRDline->j) - for (k = 0; k < atmos.Nspace; k++) Pj[k] += line->Rij[k]; - } - for (kr = 0; kr < atom->Ncont; kr++) { - continuum = &atom->continuum[kr]; - if (continuum->j == PRDline->j) - for (k = 0; k < atmos.Nspace; k++) Pj[k] += continuum->Rji[k]; - if (continuum->i == PRDline->j) - for (k = 0; k < atmos.Nspace; k++) Pj[k] += continuum->Rij[k]; - } - for (k = 0; k < atmos.Nspace; k++) { - gamma[k] = atom->n[PRDline->i][k] / atom->n[PRDline->j][k] * - PRDline->Bij / Pj[k]; + XRD = (AtomicLine **) malloc(Nsubordinate * sizeof(AtomicLine*)); + XRD[0] = PRDline; + if (input.XRD) { + for (kxrd = 0; kxrd < PRDline->Nxrd; kxrd++) + XRD[kxrd + 1] = PRDline->xrd[kxrd]; } - /* --- Store line-of-sight velocity in Doppler units to avoid - having to recompute it for every wavelength -- ------------- */ + /* --- Initialize the emission profile ratio rho -- -------------- */ - for (mu = 0; mu < atmos.Nrays; mu++) { - for (k = 0; k < atmos.Nspace; k++) - v_los[mu][k] = vproject(k, mu) / atom->vbroad[k]; + for (la = 0; la < PRDline->Nlambda; la++) { + for (mu = 0; mu < Nrays; mu++) { + for (to_obs = 0; to_obs <= 1; to_obs++) { + lamu = 2*(Nrays*la + mu) + to_obs; + for (k = 0; k < atmos.Nspace; k++) { + PRDline->rho_prd[lamu][k] = 1.0; + } + } + } } - /* --- Depth-dependent damping parameter -- -------------- */ + Pj = (double *) malloc (atmos.Nspace * sizeof(double)); + adamp = (double *) malloc (atmos.Nspace * sizeof(double)); + cDop = (NM_TO_M * PRDline->lambda0) / (4.0 * PI); - for (k = 0; k < atmos.Nspace; k++) { - Jbar[k] = PRDline->Rij[k] / PRDline->Bij; - sv[k] = 1.0 / (SQRTPI * atom->vbroad[k]); - adamp[k] = - (PRDline->Grad + PRDline->Qelast[k]) * cDop / atom->vbroad[k]; - } - /* --- Outer loop over emission wavelength -- -------------- */ + for (k = 0; k < atmos.Nspace; k++) { + adamp[k] = (PRDline->Grad + PRDline->Qelast[k]) * cDop / atom->vbroad[k]; - for (la = 0; la < PRDline->Nlambda; la++) { + /* --- Evaluate the total rate Pj out of the line's upper level - */ - /* --- Loop over emission angle (down and up) -- -------------- */ + Pj[k] = PRDline->Qelast[k]; + for (ip = 0; ip < atom->Nlevel; ip++) { + ij = ip * atom->Nlevel + PRDline->j; + Pj[k] += atom->C[ij][k]; + } + for (kr = 0; kr < atom->Nline; kr++) { + line = &atom->line[kr]; + if (line->j == PRDline->j) Pj[k] += line->Rji[k]; + if (line->i == PRDline->j) Pj[k] += line->Rij[k]; + } + for (kr = 0; kr < atom->Ncont; kr++) { + continuum = &atom->continuum[kr]; + if (continuum->j == PRDline->j) Pj[k] += continuum->Rji[k]; + if (continuum->i == PRDline->j) Pj[k] += continuum->Rij[k]; + } + } + /* --- Loop over subordinate lines -- -------------- */ - for (mu = 0; mu < atmos.Nrays; mu++) { - for (to_obs = 0; to_obs <= 1; to_obs++) { - lamu = 2*(atmos.Nrays*la + mu) + to_obs; + for (kxrd = 0; kxrd < Nsubordinate; kxrd++) { + XRDline = XRD[kxrd]; + waveratio = XRDline->lambda0 / PRDline->lambda0; - if (atmos.moving || - (PRDline->polarizable && input.StokesMode > FIELD_FREE)) - phi_emit = PRDline->phi[lamu]; - else - phi_emit = PRDline->phi[la]; + J_k = realloc(J_k, XRDline->Nlambda * sizeof(double)); + q_abs = realloc(q_abs, XRDline->Nlambda * sizeof(double)); - for (k = 0; k < atmos.Nspace; k++) { - PRDline->rho_prd[lamu][k] = 0.0; - RIInorm[k] = 0.0; - } - /* --- Depth-dependent emission wavelengths for this - direction (in Doppler units). -- -------------- */ + /* --- Loop over all spatial locations -- -------------- */ - if (to_obs) { - for (k = 0; k < atmos.Nspace; k++) - v_emit[k] = (PRDline->lambda[la] - PRDline->lambda0) * CLIGHT / - (atom->vbroad[k] * PRDline->lambda0) + v_los[mu][k]; - } else { - for (k = 0; k < atmos.Nspace; k++) - v_emit[k] = (PRDline->lambda[la] - PRDline->lambda0) * CLIGHT / - (atom->vbroad[k] * PRDline->lambda0) - v_los[mu][k]; - } - - /* --- Loop over absorption directions -- -------------- */ + for (k = 0; k < atmos.Nspace; k++) { + gamma = atom->n[XRDline->i][k] / atom->n[PRDline->j][k] * + XRDline->Bij / Pj[k]; + Jbar = XRDline->Rij[k] / XRDline->Bij; - for (mup = 0; mup < atmos.Nrays; mup++) { - wmup = 0.5 * atmos.wmu[mup]; - for (to_obs_p = 0; to_obs_p <= 1; to_obs_p++) { + /* --- Get local mean intensity and wavelength in Doppler units */ - /* --- Read specific intensity in this direction for all - wavelengths -- -------------- */ + for (la = 0; la < XRDline->Nlambda; la++) { + // Tiago: this is the line where we need to replace spectrum.J by spectrum.Jgas - for (lap = 0; lap < PRDline->Nlambda; lap++) - readImu(PRDline->Nblue + lap, mup, to_obs_p, Imup[lap]); + J_k[la] = spectrum.Jgas[XRDline->Nblue + la][k]; + //J_k[la] = spectrum.J[XRDline->Nblue + la][k]; + q_abs[la] = (XRDline->lambda[la] - XRDline->lambda0) * CLIGHT / + (XRDline->lambda0 * atom->vbroad[k]); + } + switch (representation) { + case LINEAR: + break; + case SPLINE: + splineCoef(XRDline->Nlambda, q_abs, J_k); + break; + case EXP_SPLINE: + exp_splineCoef(XRDline->Nlambda, q_abs, J_k, TENSION); + break; + } + /* --- Outer wavelength loop over emission wavelengths -- ----- */ - /* --- Loop over space -- -------------- */ + for (la = 0; la < PRDline->Nlambda; la++) { + q_emit = (PRDline->lambda[la] - PRDline->lambda0) * CLIGHT / + (PRDline->lambda0 * atom->vbroad[k]); - for (k = 0; k < atmos.Nspace; k++) { - for (lap = 0; lap < PRDline->Nlambda; lap++) { - Ik[lap] = Imup[lap][k]; + /* --- Establish integration limits over absorption wavelength, + using only regions where the redistribution function is + non-zero. (See also function GII.) -- -------------- */ + + if (fabs(q_emit) < PRD_QCORE) { + q0 = -PRD_QWING; + qN = PRD_QWING; + } else { + if (fabs(q_emit) < PRD_QWING) { + if (q_emit > 0.0) { + q0 = -PRD_QWING; + qN = waveratio * (q_emit + PRD_QSPREAD); + } else { + q0 = waveratio * (q_emit - PRD_QSPREAD); + qN = PRD_QWING; + } + } else { + q0 = waveratio * (q_emit - PRD_QSPREAD); + qN = waveratio * (q_emit + PRD_QSPREAD); + } + } + Np = (int) ((qN - q0) / PRD_DQ) + 1; + qp = (double *) realloc(qp, Np * sizeof(double)); + for (lap = 1, qp[0] = q0; lap < Np; lap++) + qp[lap] = qp[lap - 1] + PRD_DQ; - /* --- Array of absorption wavelengths in this - direction -- -------------- */ + /* --- Fold interpolation of J for symmetric lines. -- ------ */ - if (to_obs_p) - v_abs[lap] = - (PRDline->lambda[lap] - PRDline->lambda0) * CLIGHT / - (PRDline->lambda0 * atom->vbroad[k]) + v_los[mup][k]; - else - v_abs[lap] = - (PRDline->lambda[lap] - PRDline->lambda0) * CLIGHT / - (PRDline->lambda0 * atom->vbroad[k]) - v_los[mup][k]; - } - /* --- Setup spline interpolation coefficients for the - integration over absorption intensity -- ------- */ + if (XRDline->symmetric) { + qpp = (double *) realloc(qpp, Np * sizeof(double)); + for (lap = 0; lap < Np; lap++) qpp[lap] = fabs(qp[lap]); + } - switch (representation) { - case SPLINE: - splineCoef(PRDline->Nlambda, v_abs, Ik); - break; - case EXP_SPLINE: - exp_splineCoef(PRDline->Nlambda, v_abs, Ik, TENSION); - break; - case LINEAR: break; - } - /* --- Establish integration limits over absorption - wavelength, using only regions where the - redistribution function is non-zero. -- -------- */ - - if (fabs(v_emit[k]) < PRD_QCORE) { - v0 = -PRD_QWING; - vN = PRD_QWING; - } else { - if (fabs(v_emit[k]) < PRD_QWING) { - if (v_emit[k] > 0.0) { - v0 = -PRD_QWING; - vN = v_emit[k] + PRD_QSPREAD; - } else { - v0 = v_emit[k] - PRD_QSPREAD; - vN = PRD_QWING; - } - } else { - v0 = v_emit[k] - PRD_QSPREAD; - vN = v_emit[k] + PRD_QSPREAD; - } - } - Np = (int) ((vN - v0) / PRD_DQ) + 1; - vp = (double *) realloc(vp, Np * sizeof(double)); - for (lap = 1, vp[0] = v0; lap < Np; lap++) - vp[lap] = vp[lap - 1] + PRD_DQ; + /* --- Interpolate mean intensity onto fine grid. Choose linear, + spline or exponential spline interpolation -- -------- */ - /* --- Interpolate specific intensity onto fine grid. - Choose linear, spline or exponential spline - interpolation -- -------------- */ + J = (double *) realloc(J, Np * sizeof(double)); + switch (representation) { + case LINEAR: + Linear(XRDline->Nlambda, q_abs, J_k, Np, + (XRDline->symmetric) ? qpp : qp, J, hunt=TRUE); + break; + case SPLINE: + splineEval(Np, (XRDline->symmetric) ? qpp : qp, J, hunt=TRUE); + break; + case EXP_SPLINE: + exp_splineEval(Np, (XRDline->symmetric) ? qpp : qp, J, hunt=TRUE); + break; + } + /* --- Compute the redistribution weights -- -------------- */ - I = (double *) realloc(I, Np * sizeof(double)); - switch (representation) { - case LINEAR: - Linear(PRDline->Nlambda, v_abs, Ik, Np, vp, I, hunt=TRUE); - break; - case SPLINE: - splineEval(Np, vp, I, hunt=TRUE); - break; - case EXP_SPLINE: - exp_splineEval(Np, vp, I, hunt=TRUE); - break; - } - /* --- Compute the redistribution weights -- ---------- */ - - rii = (double *) realloc(rii, Np * sizeof(double)); - if (initialize) { + gii = (double *) realloc(gii, 2*Nrays*Np * sizeof(double)); + //gii = (double *) realloc(gii, Np * sizeof(double)); + + if (initialize) { - /* --- Integration weights (See: Press et al. - Numerical Recipes, p. 107, eq. 4.1.12) -- ---- */ + /* --- Integration weights (See: Press et al. Numerical Recipes, + p. 107, eq. 4.1.12) -- -------------- */ - wv = (double *) realloc(wv, Np * sizeof(double)); - wv[0] = 5.0/12.0 * PRD_DQ; - wv[1] = 13.0/12.0 * PRD_DQ; - for (lap = 2; lap < Np-2; lap++) wv[lap] = PRD_DQ; - wv[Np-1] = 5.0/12.0 * PRD_DQ; - wv[Np-2] = 13.0/12.0 * PRD_DQ; + wq = (double *) realloc(wq, Np * sizeof(double)); + wq[0] = 5.0/12.0 * PRD_DQ; + wq[1] = 13.0/12.0 * PRD_DQ; + for (lap = 2; lap < Np-2; lap++) wq[lap] = PRD_DQ; + wq[Np-1] = 5.0/12.0 * PRD_DQ; + wq[Np-2] = 13.0/12.0 * PRD_DQ; - /* --- RII(x, mu, x', mu')/phi(x, mu) -- ------------ */ + // Tiago: adding angle loops here + for (mu = 0; mu < Nrays; mu++) { + for (to_obs = 0; to_obs <= 1; to_obs++) { + + sign = (to_obs) ? 1.0 : -1.0; + + //vlos = sign*spectrum.v_los[mu][k] / atom->vbroad[k]; + vlos = 0.; + + for (lap = 0; lap < Np; lap++) { + + idx = mu*2*Np + to_obs*Np + lap; /* index of gii */ - for (lap = 0; lap < Np; lap++) { - rii[lap] = RII(v_emit[k], vp[lap], adamp[k], mu, mup) * - (sv[k] / phi_emit[k]) * wv[lap] * wmup; - } - if ((Nwrite = fwrite(rii, sizeof(double), Np, - PRDline->fp_GII)) != Np) { - sprintf(messageStr, + gii[idx] = GII(adamp[k], waveratio, q_emit-vlos, qp[lap]) * wq[lap]; + } + } + } + + if ((Nwrite = + //fwrite(gii, sizeof(double), Np, PRDline->fp_GII)) != Np) { + fwrite(gii, sizeof(double), 2*Nrays*Np, PRDline->fp_GII)) != 2*Nrays*Np) { + sprintf(messageStr, "Unable to write proper number of redistribution weights\n" " Wrote %d instead of %d.\n Line %d -> %d, la = %d, k = %d", - Nwrite, Np, PRDline->j, PRDline->i, la, k); - Error(ERROR_LEVEL_2, routineName, messageStr); - } - } else { - if ((Nread = fread(rii, sizeof(double), Np, - PRDline->fp_GII)) != Np) { - sprintf(messageStr, - "Unable to read proper number of redistribution weights\n" - " Read %d instead of %d.\n Line %d -> %d, la = %d, k = %d", - Nread, Np, PRDline->j, PRDline->i, la, k); - Error(ERROR_LEVEL_2, routineName, messageStr); - } - } - /* --- Inner wavelength loop doing actual wavelength - integration over absorption wavelengths -- ----- */ + Nwrite, Np, XRDline->j, XRDline->i, la, k); + Error(ERROR_LEVEL_2, routineName, messageStr); + } + } else { + if ((Nread = + //fread(gii, sizeof(double), Np, PRDline->fp_GII)) != Np) { + fread(gii, sizeof(double), 2*Nrays*Np, PRDline->fp_GII)) != 2*Nrays*Np) { + sprintf(messageStr, + "Unable to read proper number of redistribution weights\n" + " Read %d instead of %d.\n Line %d -> %d, la = %d, k = %d", + Nread, Np, XRDline->j, XRDline->i, la, k); + Error(ERROR_LEVEL_2, routineName, messageStr); + } + } + /* --- Inner wavelength loop doing actual wavelength integration + over absorption wavelengths -- -------------- */ - for (lap = 0; lap < Np; lap++) { - RIInorm[k] += rii[lap]; - PRDline->rho_prd[lamu][k] += I[lap] * rii[lap]; - } // inner wavelength loop - } // depth loop - } // direction loop - } // mup loop + for (mu = 0; mu < Nrays; mu++) { + for (to_obs = 0; to_obs <= 1; to_obs++) { + + lamu = 2*(Nrays*la + mu) + to_obs; + //lamu = la; + + gnorm = 0.0; + scatInt = 0.0; + + // Integral in absorption wavelength: + for (lap = 0; lap < Np; lap++) { + idx = mu*2*Np + to_obs*Np + lap; /* index of gii */ + //idx = lap; + gnorm += gii[idx]; + scatInt += J[lap] * gii[idx]; + } + + // Tiago: put condition to prevent negative scatInt? + PRDline->rho_prd[lamu][k] += gamma*(scatInt/gnorm - Jbar); /* Integral over XRD lines*/ + } /* end of to_obs loop */ + } /* end of mu loop */ - for (k = 0; k < atmos.Nspace; k++) { - PRDline->rho_prd[lamu][k] = 1.0 + - gamma[k] * (PRDline->rho_prd[lamu][k]/RIInorm[k] - Jbar[k]); + // Tiago: must interpolate rho_prd back into the lab frame. use Linear() + } /* End of emission wavelength loop */ + } /* End of depth loop */ + } /* End of subordinate line loop */ + + /* --- Tiago: shift rho back to lab rest frame -- -------------- */ + + rho_stor = (double *) calloc(PRDline->Nlambda, sizeof(double)); + rho_intp = (double *) calloc(PRDline->Nlambda, sizeof(double)); + lambda_s = (double *) calloc(PRDline->Nlambda, sizeof(double)); + rho_new = matrix_double(2*atmos.Nrays * PRDline->Nlambda, atmos.Nspace); + + /* + for (k = 0; k < atmos.Nspace; k++) { + for (mu = 0; mu < Nrays; mu++ ) { + for (to_obs = 0; to_obs <= 1; to_obs++) { + + // Grab rho[lambda] and shifted wavelengths + for (la = 0; la < PRDline->Nlambda; la++) { + lamu = 2*(Nrays*la + mu) + to_obs; + sign = (to_obs) ? 1.0 : -1.0; + rho_stor[la] = PRDline->rho_prd[lamu][k]; + lambda_s[la] = PRDline->lambda[la]*(1.-spectrum.v_los[mu][k]*sign/CLIGHT); + } + + // Linearly interpolate, equivalent to wavelength shift + Linear(PRDline->Nlambda, lambda_s, rho_stor, PRDline->Nlambda, + PRDline->lambda, rho_intp, hunt=TRUE); + + // Save back into rho + for (la = 0; la < PRDline->Nlambda; la++) { + lamu = 2*(Nrays*la + mu) + to_obs; + rho_new[lamu][k] = rho_intp[la]; + PRDline->rho_prd[lamu][k] = rho_intp[la]; } - } - } - } + } // direction loop + } // mu loop + + } // depth loop + */ + + + free(rho_stor); free(rho_intp); free(lambda_s); freeMatrix((void **) rho_new); + /* --- Clean temporary variable space -- -------------- */ - freeMatrix((void **) Imup); - freeMatrix((void **) v_los); - - free(Ik); free(v_abs); free(v_emit); free(adamp); - free(gamma); free(wv); free(I); free(rii); - free(vp); free(Jbar); free(RIInorm); free(sv); + for (kxrd = 0; kxrd < Nsubordinate; kxrd++) { + if (XRD[kxrd]->symmetric) { + free(qpp); + break; + } + } + free(J_k); free(q_abs); free(gii); + free(qp); free(wq); free(J); + free(XRD); free(Pj); free(adamp); sprintf(messageStr, "Scatter Int %5.1f", PRDline->lambda0); getCPU(3, TIME_POLL, messageStr); } -/* ------- end ---------------------------- PRDAngleScatter.c ------- */ +/* ------- end ---------------------------- PRDAngleApproxScatter.c - */