94 long int ion , nshell , low , ihi , ipop , ip;
101 if( !
dense.lgElmtOn[nelem] )
111 nshell =
Heavy.nsShells[nelem][ion] - 1;
114 low =
opac.ipElement[nelem][ion][nshell][0];
115 ihi =
opac.ipElement[nelem][ion][nshell][1];
116 ipop =
opac.ipElement[nelem][ion][nshell][2];
119 for( ip=low-1; ip < ihi; ip++ )
121 opac.OpacStack[ip-low+ipop] *= scale;
147 static const realnum csh2p[
NCSH2P]={6.75f,0.24f,8.68f,2.5f,10.54f,7.1f,12.46f,
167 fprintf(
ioQQQ,
" OpacityCreateAll called but NOT evaluated since already done.\n" );
178 fprintf(
ioQQQ,
" OpacityCreateAll called, evaluating.\n" );
182 for( i=0; i <
rfield.nupper; i++ )
184 opac.opacity_abs[i] = 0.;
193 for( nelem=ipISO; nelem <
LIMELM; nelem++ )
195 iso_sp[ipISO][nelem].HighestLevelOpacStack.resize(0);
196 if(
dense.lgElmtOn[nelem] )
203 opac.ipElement[nelem][nelem-ipISO][0][2] =
opac.nOpacTot + 1;
210 for(
long index=0; index <
iso_sp[ipISO][nelem].numLevels_max; index++ )
213 iso_sp[ipISO][nelem].fb[index].ipOpac =
opac.nOpacTot + 1;
219 iso_sp[ipISO][nelem].fb[index].xIsoLevNIonRyd );
222 need = nupper -
iso_sp[ipISO][nelem].fb[index].ipIsoLevNIonCon + 1;
228 for( i=
iso_sp[ipISO][nelem].fb[index].ipIsoLevNIonCon-1; i < nupper; i++ )
232 opac.OpacStack[i-
iso_sp[ipISO][nelem].fb[index].ipIsoLevNIonCon+
iso_sp[ipISO][nelem].fb[index].ipOpac] = crs;
233 if( index==
iso_sp[ipISO][nelem].numLevels_max-1 )
234 iso_sp[ipISO][nelem].HighestLevelOpacStack.push_back( crs );
237 opac.nOpacTot += need;
246 if( (*diatom)->lgEnabled &&
mole_global.lgStancil )
248 for( vector< diss_tran >::iterator tran = (*diatom)->Diss_Trans.begin(); tran != (*diatom)->Diss_Trans.end(); ++tran )
251 long lower_limit =
ipoint(tran->energies[0]);
252 long upper_limit =
ipoint(tran->energies.back());
253 upper_limit =
MIN2( upper_limit,
rfield.nflux-1 );
254 long ipCont_Diss =
opac.nOpacTot + 1;
261 for(i = lower_limit; i <= upper_limit; ++i)
263 opac.OpacStack[ipCont_Diss + num_points - 1] =
267 opac.nOpacTot += num_points;
291 for( i=0; i <
opac.ipCKshell; i++ )
293 opac.OpacStack[i-1+
opac.iopcom] = thom;
300 for( i=
opac.ipCKshell; i <
rfield.nupper; i++ )
302 dx =
rfield.AnuOrg[i]/3.7573e4;
304 opac.OpacStack[i-1+
opac.iopcom] = thom*3.e0/4.e0*((1.e0 +
305 dx)/
POW3(dx)*(2.e0*dx*(1.e0 + dx)/(1.e0 + 2.e0*dx) - log(1.e0+
306 2.e0*dx)) + 1.e0/2.e0/dx*log(1.e0+2.e0*dx) - (1.e0 + 3.e0*
307 dx)/
POW3(1.e0 + 2.e0*dx));
328 x =
rfield.AnuOrg[i]/7.512e4*2.;
331 POW2((-0.46737118 + x*(0.349255416 + x*0.002179893))/(1. +
332 x*(0.130471301 + x*0.000524906)));
341 for( i=0; i <
rfield.nupper; i++ )
354 for( i=0; i <
rfield.nupper; i++ )
379 (*diatom)->ip_photo_opac_offset =
opac.nOpacTot + 1;
380 opac.nOpacTot += (*diatom)->OpacityCreate(
opac.OpacStack );
410 for( nelem=2; nelem <
LIMELM; nelem++ )
412 if(
dense.lgElmtOn[nelem] )
445 crs =
ofit(eps,opart);
449 for(
long n=1; n < 6; n++ )
485 crs = 3.85e-18*(4.4*pow(
rfield.AnuOrg[i]/thres,-1.5) - 3.38*
486 pow(
rfield.AnuOrg[i]/thres,-2.5));
499 (0.2602325880970085 +
500 445.8558249365131*exp(-
rfield.AnuOrg[i]/0.1009243952792674))*
514 " OpacityCreateAll return OK, number of opacity cells used in OPSC= %ld and OPSV is dimensioned %ld\n",
520 if(
opac.lgCompileOpac )
522 fprintf(
ioQQQ,
"The COMPILE OPACITIES command is currently not supported\n" );
528 " Please consider revising ndimOpacityStack in OpacityCreateAll, a total of %li cells were needed.\n\n" ,
opac.nOpacTot);
553 *ip =
opac.nOpacTot + 1;
556 thres =
rfield.anu[ilo-1];
561 for( i=ilo-1; i < ihi; i++ )
563 opac.OpacStack[i-ilo+*ip] = cross*pow(
rfield.anu[i]/thres,-s);
566 opac.nOpacTot += ihi - ilo + 1;
594 *ipop =
opac.nOpacTot + 1;
600 fprintf(
ioQQQ,
" Too many opacities were entered into OpacityCreateReilMan. Increase the value of NOP.\n" );
601 fprintf(
ioQQQ,
" chLabl was %4.4s\n", chLabl );
608 for( i=0; i < ncr; i++ )
610 en[i] = cross[i*2]/13.6f;
611 cs[i] = cross[(i+1)*2-1]*1e-18f;
615 if( en[0] >
rfield.anu[low-1] )
618 " OpacityCreateReilMan: The entered opacity energy bandwidth is not large enough (low fail).\n" );
620 " The desired energy (Ryd) was%12.5eeV and the lowest entered in the array was%12.5e eV\n",
623 fprintf(
ioQQQ,
" chLabl was %4.4s\n", chLabl );
624 fprintf(
ioQQQ,
" The original energy (eV) and cross section (mb) arrays follow:\n" );
625 fprintf(
ioQQQ,
" " );
627 for( i=0; i < ncross; i++ )
629 fprintf(
ioQQQ,
"%11.4e", cross[i] );
632 fprintf(
ioQQQ,
"\n" );
636 slope = (cs[1] - cs[0])/(en[1] - en[0]);
643 for( i=low-1; i < ihi; i++ )
645 if(
rfield.anu[i] > en[ics-1] &&
rfield.anu[i] <= en[ics] )
647 opac.OpacStack[i-low+*ipop] = cs[ics-1] + slope*(
rfield.anu[i] -
656 fprintf(
ioQQQ,
" OpacityCreateReilMan: The entered opacity energy bandwidth is not large enough (high fail).\n" );
657 fprintf(
ioQQQ,
" The entered energy was %10.2eeV and the highest in the array was %10.2eeV\n",
658 rfield.anu[i]*13.6, en[ncr-1]*13.6 );
659 fprintf(
ioQQQ,
" chLabl was %4.4s\n", chLabl
661 fprintf(
ioQQQ,
" The lowest energy enterd in the array was%10.2e eV\n",
663 fprintf(
ioQQQ,
" The highest energy ever needed would be%10.2eeV\n",
665 fprintf(
ioQQQ,
" The lowest energy needed was%10.2eeV\n",
670 slope = (cs[ics] - cs[ics-1])/(en[ics] - en[ics-1]);
671 if(
rfield.anu[i] > en[ics-1] &&
rfield.anu[i] <= en[ics] )
673 opac.OpacStack[i-low+*ipop] = cs[ics-1] + slope*(
rfield.anu[i] -
679 fprintf(
ioQQQ,
" Internal logical error in OpacityCreateReilMan.\n" );
680 fprintf(
ioQQQ,
" The desired energy (%10.2eeV), I=%5ld, is not within the next energy bound%10.2e%10.2e\n",
681 rfield.anu[i]*13.6, i, en[ics-1], en[ics] );
683 fprintf(
ioQQQ,
" The previous energy (eV) was%10.2e\n",
686 fprintf(
ioQQQ,
" Here comes the energy array. ICS=%4ld\n",
689 for( j=0; j < ncr; j++ )
691 fprintf(
ioQQQ,
"%10.2e", en[j] );
693 fprintf(
ioQQQ,
"\n" );
695 fprintf(
ioQQQ,
" chLabl was %4.4s\n", chLabl );
702 opac.nOpacTot += ihi - low + 1;
715 static const double y[7][5] = {
716 {8.915,3995.,3.242,10.44,0.0},
717 {11.31,1498.,5.27,7.319,0.0},
718 {10.5,1.059e05,1.263,13.04,0.0},
719 {19.49,48.47,8.806,5.983,0.0},
720 {50.,4.244e04,0.1913,7.012,4.454e-02},
721 {110.5,0.1588,148.3,-3.38,3.589e-02},
722 {177.4,32.37,381.2,1.083,0.0}
724 static const double eth[7]={13.62,16.94,18.79,28.48,50.,110.5,538.};
725 static const long l[7]={1,1,1,0,1,1,0};
742 for(
int i=0; i < 7; i++ )
747 for(
int i=0; i < 7; i++ )
753 q = 5.5 - 0.5*y[i][3] + l[i];
758 pow(x,q)/pow(1.0 + sqrt(x/y[i][2]),y[i][3]));
793 for( ion=0; ion < nelem; ion++ )
797 for( ip=0; ip <
rfield.nupper; ip++ )
799 opac.opacity_abs[ip] = 0.;
803 nelec = nelem+1 - ion;
806 for( nshell=0; nshell <
Heavy.nsShells[nelem][ion]; nshell++ )
809 opac.ipElement[nelem][ion][nshell][2] =
opac.nOpacTot + 1;
812 limit =
opac.ipElement[nelem][ion][nshell][1];
815 need = limit -
opac.ipElement[nelem][ion][nshell][0] + 1;
822 low =
opac.ipElement[nelem][ion][nshell][0];
823 ihi =
opac.ipElement[nelem][ion][nshell][1];
824 ipop =
opac.ipElement[nelem][ion][nshell][2];
828 ASSERT( low <= ihi || low<5 );
831 for( ip=low-1; ip < ihi; ip++ )
842 opac.OpacStack[ip-low+ipop] = cs*1e-18;
845 opac.opacity_abs[ip] += cs;
848 opac.nOpacTot += ihi - low + 1;
851 if(
save.lgPunPoint )
853 fprintf(
save.ipPoint,
"%3ld%3ld%3ld%10.2e%10.2e%10.2e%10.2e\n",
854 nelem, ion, nshell,
rfield.anu[low-1],
rfield.anu[ihi-1],
855 opac.OpacStack[ipop-1],
opac.OpacStack[ihi-low+ipop-1] );
861 for( ip=
opac.ipElement[nelem][ion][
Heavy.nsShells[nelem][ion]-1][0]-1;
880 fprintf(
ioQQQ,
" NOTE OpacityCreate1Element needed more opacity cells than ndimOpacityStack, please consider increasing it.\n" );
881 fprintf(
ioQQQ,
" NOTE OpacityCreate1Element doubled memory allocation to %li.\n",
ndimOpacityStack );
910 ASSERT( crs > 0. && crs < 1e-10 );
912 else if( index <
iso_sp[ipISO][nelem].numLevels_max -
iso_sp[ipISO][nelem].nCollapsed_max )
915 double photon =
MAX2( EgammaRyd/
iso_sp[ipISO][nelem].fb[index].xIsoLevNIonRyd, 1. + FLT_EPSILON*2. );
919 ASSERT( crs > 0. && crs < 1e-10 );
929 EgammaRyd =
MAX2( EgammaRyd ,
iso_sp[ipISO][nelem].fb[index].xIsoLevNIonRyd*1.001f );
933 ASSERT( crs > 0. && crs < 1e-10 );
938 double photon =
MAX2( EgammaRyd/
iso_sp[ipISO][nelem].fb[index].xIsoLevNIonRyd, 1. + FLT_EPSILON*2. );
946 ASSERT( crs > 0. && crs < 1e-10 );
951 EgammaRyd =
MAX2( EgammaRyd ,
iso_sp[ipISO][nelem].fb[index].xIsoLevNIonRyd);
955 long int nup =
iso_sp[
ipHE_LIKE][nelem].n_HighestResolved_max + index + 1 -
963 (EgammaRyd <
iso_sp[ipISO][nelem].fb[index].xIsoLevNIonRyd*1.02) ||
964 (crs > 0. && crs < 1e-10) );
979 ASSERT( crs > 0. && crs < 1e-10 );
996 static double y2[
NCRS];
997 static double crs[
NCRS]={0.,0.124,0.398,0.708,1.054,1.437,1.805,
998 2.176,2.518,2.842,3.126,3.377,3.580,3.741,3.851,3.913,3.925,
999 3.887,3.805,3.676,3.511,3.306,3.071,2.810,2.523,2.219,1.898,
1000 1.567,1.233,.912,.629,.39,.19};
1001 static double ener[
NCRS]={0.,0.001459,0.003296,0.005256,0.007351,
1002 0.009595,0.01201,0.01460,0.01741,0.02044,0.02375,0.02735,0.03129,
1003 0.03563,0.04043,0.04576,0.05171,0.05841,0.06601,0.07469,0.08470,
1004 0.09638,0.1102,0.1268,0.1470,0.1723,0.2049,0.2483,0.3090,0.4001,
1005 0.5520,0.8557,1.7669};
1023 energy = freq - 0.05552;
1024 if( energy < ener[0] || energy > ener[
NCRS-1] )
1057 rayleh_v = (8.41e-25*
powi(ener,4) + 3.37e-24*
powi(ener,6))*
1061 else if( ener < 0.646 )
1063 rayleh_v = (8.41e-25*
powi(ener,4) + 3.37e-24*
powi(ener,6) +
1064 4.71e-22*
powi(ener,14))*
hydro.DampOnFac;
1067 else if( ener >= 0.646 && ener < 1.0 )
1069 rayleh_v = fabs(0.74959-ener);
1072 rayleh_v =
MAX2(rayleh_v,1e-24)*
hydro.DampOnFac;
const int NHYDRO_MAX_LEVEL
NORETURN void TotalInsanity(void)
#define DEBUG_ENTRY(funcname)
double powi(double, long int)
double phfit(long int nz, long int ne, long int is, double e)
double hpfit(long int iz, long int n, double e)
long ipoint(double energy_ryd)
vector< diatomics * > diatoms
vector< diatomics * >::iterator diatom_iter
double MolDissocCrossSection(const diss_tran &tran, const double &Mol_Ene)
double He_cross_section(double EgammaRyd, double EthRyd, long n, long l, long S, long nelem)
double H_photo_cs(double rel_photon_energy, long int n, long int l, long int iz)
t_iso_sp iso_sp[NISO][LIMELM]
t_mole_global mole_global
STATIC void opacity_more_memory(void)
STATIC void OpacityCreatePowerLaw(long int ilo, long int ihi, double cross, double s, long int *ip)
STATIC double rayleh(double ener)
STATIC void OpacityCreate1Element(long int nelem)
STATIC double Opacity_iso_photo_cs(double energy, long ipISO, long nelem, long index)
STATIC double hmiopc(double freq)
STATIC void OpacityCreateReilMan(long int low, long int ihi, const realnum cross[], long int ncross, long int *ipop, const char *chLabl)
STATIC double ofit(double e, realnum opart[])
static long int ndimOpacityStack
STATIC void OpacityValenceRescale(long int nelem, double scale)
void OpacityCreateAll(void)
UNUSED const double FR1RYD
UNUSED const double EVRYD
void spline(const double x[], const double y[], long int n, double yp1, double ypn, double y2a[])
void splint(const double xa[], const double ya[], const double y2a[], long int n, double x, double *y)