00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #include "cddefines.h"
00016 #include "physconst.h"
00017 #include "dense.h"
00018 #include "continuum.h"
00019 #include "iso.h"
00020 #include "hydrogenic.h"
00021 #include "oxy.h"
00022 #include "trace.h"
00023 #include "heavy.h"
00024 #include "rfield.h"
00025 #include "hmi.h"
00026 #include "atmdat.h"
00027 #include "save.h"
00028 #include "grains.h"
00029 #include "thirdparty.h"
00030 #include "hydro_bauman.h"
00031 #include "opacity.h"
00032 #include "helike_recom.h"
00033 #include "taulines.h"
00034 #include "h2.h"
00035 #include "h2_priv.h"
00036 #include "ipoint.h"
00037 #include "mole.h"
00038
00039 static const int NCSH2P = 10;
00040
00041
00042 static long int ndimOpacityStack = 2600000L;
00043
00044
00045 STATIC void OpacityCreate1Element(long int nelem);
00046
00047
00048 STATIC void opacity_more_memory(void);
00049
00050
00051 STATIC double hmiopc(double freq);
00052
00053
00054 STATIC double rayleh(double ener);
00055
00056
00057 STATIC double Opacity_iso_photo_cs( double energy , long ipISO , long nelem , long index );
00058
00059
00060 STATIC void OpacityCreateReilMan(long int low,
00061 long int ihi,
00062 const realnum cross[],
00063 long int ncross,
00064 long int *ipop,
00065 const char *chLabl);
00066
00067 static bool lgRealloc = false;
00068
00069
00070 STATIC void OpacityCreatePowerLaw(
00071
00072 long int ilo,
00073
00074 long int ihi,
00075
00076 double cross,
00077
00078 double s,
00079
00080 long int *ip);
00081
00082
00083 STATIC double ofit(double e,
00084 realnum opart[]);
00085
00086
00087 STATIC void OpacityValenceRescale(
00088
00089 long int nelem ,
00090
00091 double scale )
00092 {
00093
00094 long int ion , nshell , low , ihi , ipop , ip;
00095
00096 DEBUG_ENTRY( "OpacityValenceRescale()" );
00097
00098
00099
00100
00101 if( !dense.lgElmtOn[nelem] )
00102 {
00103 return;
00104 }
00105
00106
00107 ASSERT( scale >= 0. );
00108
00109 ion = 0;
00110
00111 nshell = Heavy.nsShells[nelem][ion] - 1;
00112
00113
00114 low = opac.ipElement[nelem][ion][nshell][0];
00115 ihi = opac.ipElement[nelem][ion][nshell][1];
00116 ipop = opac.ipElement[nelem][ion][nshell][2];
00117
00118
00119 for( ip=low-1; ip < ihi; ip++ )
00120 {
00121 opac.OpacStack[ip-low+ipop] *= scale;
00122 }
00123 return;
00124 }
00125
00126 void OpacityCreateAll(void)
00127 {
00128 long int i,
00129 ipISO ,
00130 need ,
00131 nelem;
00132
00133 realnum opart[7];
00134
00135 double crs,
00136 dx,
00137 eps,
00138 thom,
00139 thres,
00140 x;
00141
00142
00143
00144
00145
00146
00147 static const realnum csh2p[NCSH2P]={6.75f,0.24f,8.68f,2.5f,10.54f,7.1f,12.46f,
00148 6.0f,14.28f,2.7f};
00149
00150 DEBUG_ENTRY( "OpacityCreateAll()" );
00151
00152
00153
00154
00155
00156
00157 GrainsInit();
00158
00159
00160
00161
00162 if( lgOpacMalloced )
00163 {
00164
00165 if( trace.lgTrace )
00166 {
00167 fprintf( ioQQQ, " OpacityCreateAll called but NOT evaluated since already done.\n" );
00168 }
00169 return;
00170 }
00171
00172 lgOpacMalloced = true;
00173
00174
00175 opac.OpacStack = (double*)MALLOC((size_t)ndimOpacityStack*sizeof(double));
00176
00177 if( trace.lgTrace )
00178 {
00179 fprintf( ioQQQ, " OpacityCreateAll called, evaluating.\n" );
00180 }
00181
00182
00183 for( i=0; i < rfield.nupper; i++ )
00184 {
00185 opac.opacity_abs[i] = 0.;
00186 }
00187
00188
00189 opac.nOpacTot = 0;
00190
00191
00192 for( ipISO=ipH_LIKE; ipISO<=ipHE_LIKE; ++ipISO )
00193 {
00194 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00195 {
00196 if( dense.lgElmtOn[nelem] )
00197 {
00198 long int nupper;
00199
00200
00201
00202
00203 opac.ipElement[nelem][nelem-ipISO][0][2] = opac.nOpacTot + 1;
00204
00205 fixit();
00206
00207
00208
00209 nupper = rfield.nupper;
00210 for( long index=0; index < iso_sp[ipISO][nelem].numLevels_max; index++ )
00211 {
00212
00213 iso_sp[ipISO][nelem].fb[index].ipOpac = opac.nOpacTot + 1;
00214
00215
00216
00217
00218 ASSERT( rfield.AnuOrg[iso_sp[ipISO][nelem].fb[index].ipIsoLevNIonCon-1] > 0.94f *
00219 iso_sp[ipISO][nelem].fb[index].xIsoLevNIonRyd );
00220
00221
00222 need = nupper - iso_sp[ipISO][nelem].fb[index].ipIsoLevNIonCon + 1;
00223 ASSERT( need > 0 );
00224
00225 if( opac.nOpacTot + need > ndimOpacityStack )
00226 opacity_more_memory();
00227
00228 for( i=iso_sp[ipISO][nelem].fb[index].ipIsoLevNIonCon-1; i < nupper; i++ )
00229 {
00230 double crs =
00231 Opacity_iso_photo_cs( rfield.AnuOrg[i] , ipISO , nelem , index );
00232 opac.OpacStack[i-iso_sp[ipISO][nelem].fb[index].ipIsoLevNIonCon+iso_sp[ipISO][nelem].fb[index].ipOpac] = crs;
00233 if( index==iso_sp[ipISO][nelem].numLevels_max-1 )
00234 iso_sp[ipISO][nelem].HighestLevelOpacStack.push_back( crs );
00235 }
00236
00237 opac.nOpacTot += need;
00238 }
00239 }
00240 }
00241 }
00242
00243
00244 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
00245 {
00246 if( (*diatom)->lgEnabled && mole_global.lgStancil )
00247 {
00248 for( vector< diss_tran >::iterator tran = (*diatom)->Diss_Trans.begin(); tran != (*diatom)->Diss_Trans.end(); ++tran )
00249 {
00250
00251 long lower_limit = ipoint(tran->energies[0]);
00252 long upper_limit = ipoint(tran->energies.back());
00253 upper_limit = MIN2( upper_limit, rfield.nflux-1 );
00254 long ipCont_Diss = opac.nOpacTot + 1;
00255 long num_points = 0;
00256
00257
00258 if( opac.nOpacTot + upper_limit - lower_limit + 1 > ndimOpacityStack )
00259 opacity_more_memory();
00260
00261 for(i = lower_limit; i <= upper_limit; ++i)
00262 {
00263 opac.OpacStack[ipCont_Diss + num_points - 1] =
00264 MolDissocCrossSection( *tran, rfield.anu[i] );
00265 ++num_points;
00266 }
00267 opac.nOpacTot += num_points;
00268 }
00269 }
00270 }
00271
00272
00273 if( opac.nOpacTot + iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon + rfield.nupper > ndimOpacityStack )
00274 opacity_more_memory();
00275
00276
00277 opac.ipRayScat = opac.nOpacTot + 1;
00278 for( i=0; i < iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon; i++ )
00279 {
00280 opac.OpacStack[i-1+opac.ipRayScat] = rayleh(rfield.AnuOrg[i]);
00281 }
00282 opac.nOpacTot += iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon;
00283
00284
00285
00286
00287
00288
00289 thom = 6.65e-25;
00290 opac.iopcom = opac.nOpacTot + 1;
00291 for( i=0; i < opac.ipCKshell; i++ )
00292 {
00293 opac.OpacStack[i-1+opac.iopcom] = thom;
00294
00295
00296 }
00297
00298
00299
00300 for( i=opac.ipCKshell; i < rfield.nupper; i++ )
00301 {
00302 dx = rfield.AnuOrg[i]/3.7573e4;
00303
00304 opac.OpacStack[i-1+opac.iopcom] = thom*3.e0/4.e0*((1.e0 +
00305 dx)/POW3(dx)*(2.e0*dx*(1.e0 + dx)/(1.e0 + 2.e0*dx) - log(1.e0+
00306 2.e0*dx)) + 1.e0/2.e0/dx*log(1.e0+2.e0*dx) - (1.e0 + 3.e0*
00307 dx)/POW3(1.e0 + 2.e0*dx));
00308
00309
00310 }
00311 opac.nOpacTot += rfield.nupper - 1 + 1;
00312
00313
00314
00315
00316
00317 if( opac.nOpacTot + 3*rfield.nupper - opac.ippr + iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon - hmi.iphmin + 2 > ndimOpacityStack )
00318 opacity_more_memory();
00319
00320
00321 opac.ioppr = opac.nOpacTot + 1;
00322 for( i=opac.ippr-1; i < rfield.nupper; i++ )
00323 {
00324
00325
00326
00327
00328 x = rfield.AnuOrg[i]/7.512e4*2.;
00329
00330 opac.OpacStack[i-opac.ippr+opac.ioppr] = 5.793e-28*
00331 POW2((-0.46737118 + x*(0.349255416 + x*0.002179893))/(1. +
00332 x*(0.130471301 + x*0.000524906)));
00333
00334
00335 }
00336 opac.nOpacTot += rfield.nupper - opac.ippr + 1;
00337
00338
00339 opac.ipBrems = opac.nOpacTot + 1;
00340
00341 for( i=0; i < rfield.nupper; i++ )
00342 {
00343
00344
00345 opac.OpacStack[i-1+opac.ipBrems] =
00346
00347
00348
00349 1.03680e-8/POW3(rfield.AnuOrg[i]);
00350 }
00351 opac.nOpacTot += rfield.nupper - 1 + 1;
00352
00353 opac.iphmra = opac.nOpacTot + 1;
00354 for( i=0; i < rfield.nupper; i++ )
00355 {
00356
00357 opac.OpacStack[i-1+opac.iphmra] = 0.1175*rfield.anusqr[i];
00358 }
00359 opac.nOpacTot += rfield.nupper - 1 + 1;
00360
00361 opac.iphmop = opac.nOpacTot + 1;
00362 for( i=hmi.iphmin-1; i < iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon; i++ )
00363 {
00364
00365 opac.OpacStack[i-hmi.iphmin+opac.iphmop] =
00366 hmiopc(rfield.AnuOrg[i]);
00367 }
00368 opac.nOpacTot += iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon - hmi.iphmin + 1;
00369
00370
00371
00372
00373
00374 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
00375 {
00376 if( opac.nOpacTot + rfield.nupper - (*diatom)->ip_photo_opac_thresh + 1 > ndimOpacityStack )
00377 opacity_more_memory();
00378
00379 (*diatom)->ip_photo_opac_offset = opac.nOpacTot + 1;
00380 opac.nOpacTot += (*diatom)->OpacityCreate( opac.OpacStack );
00381 }
00382
00383
00384
00385
00386 OpacityCreateReilMan(opac.ih2pnt[0],opac.ih2pnt[1],csh2p,NCSH2P,&opac.ih2pof,
00387 "H2+ ");
00388
00389
00390
00391 if( opac.nOpacTot + rfield.nupper - iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon + 1 > ndimOpacityStack )
00392 opacity_more_memory();
00393
00394
00395 opac.iophe1[0] = opac.nOpacTot + 1;
00396 opac.ipElement[ipHELIUM][0][0][2] = opac.iophe1[0];
00397 for( i=iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon-1; i < rfield.nupper; i++ )
00398 {
00399 crs = t_ADfA::Inst().phfit(2,2,1,rfield.AnuOrg[i]*EVRYD);
00400 opac.OpacStack[i-iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon+opac.iophe1[0]] =
00401 crs*1e-18;
00402 }
00403 opac.nOpacTot += rfield.nupper - iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon + 1;
00404
00405
00406
00407
00408
00409
00410 for( nelem=2; nelem < LIMELM; nelem++ )
00411 {
00412 if( dense.lgElmtOn[nelem] )
00413 {
00414 OpacityCreate1Element(nelem);
00415 }
00416 }
00417
00418
00419
00420
00421
00422
00423 OpacityValenceRescale( ipPOTASSIUM , 5. );
00424
00425
00426
00427
00428
00429 OpacityCreatePowerLaw(opac.in1[0],opac.in1[1],9e-18,1.75,&opac.in1[2]);
00430
00431
00432
00433 if( dense.lgElmtOn[ipOXYGEN] && t_ADfA::Inst().get_version() == PHFIT96 )
00434 {
00435
00436
00437 if( opac.nOpacTot + opac.ipElement[ipOXYGEN][0][2][1] - opac.ipElement[ipOXYGEN][0][2][0] + 1 > ndimOpacityStack )
00438 opacity_more_memory();
00439
00440
00441 for( i=opac.ipElement[ipOXYGEN][0][2][0]-1; i < opac.ipElement[ipOXYGEN][0][2][1]; i++ )
00442 {
00443
00444 eps = rfield.AnuOrg[i]*EVRYD;
00445 crs = ofit(eps,opart);
00446
00447
00448 crs = opart[0];
00449 for( long n=1; n < 6; n++ )
00450 {
00451
00452 crs += opart[n];
00453 }
00454
00455 crs *= 1e-18;
00456 opac.OpacStack[i-opac.ipElement[ipOXYGEN][0][2][0]+opac.ipElement[ipOXYGEN][0][2][2]] = crs;
00457 }
00458
00459
00460 opac.nOpacTot += opac.ipElement[ipOXYGEN][0][2][1] - opac.ipElement[ipOXYGEN][0][2][0] + 1;
00461 }
00462
00463
00464 OpacityCreatePowerLaw(opac.ipo1exc[0],opac.ipo1exc[1],4.64e-18,0.,&opac.ipo1exc[2]);
00465
00466
00467
00468 OpacityCreatePowerLaw(opac.ipo3exc[0],opac.ipo3exc[1],3.8e-18,0.,&opac.ipo3exc[2]);
00469
00470
00471 OpacityCreatePowerLaw(opac.ipo3exc3[0],opac.ipo3exc3[1],5.5e-18,0.01,
00472 &opac.ipo3exc3[2]);
00473
00474
00475
00476 if( opac.nOpacTot + iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon - oxy.i2d + 1
00477 + iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon - opac.ipmgex + 1 > ndimOpacityStack )
00478 opacity_more_memory();
00479
00480
00481 opac.iopo2d = opac.nOpacTot + 1;
00482 thres = rfield.AnuOrg[oxy.i2d-1];
00483 for( i=oxy.i2d-1; i < iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon; i++ )
00484 {
00485 crs = 3.85e-18*(4.4*pow(rfield.AnuOrg[i]/thres,-1.5) - 3.38*
00486 pow(rfield.AnuOrg[i]/thres,-2.5));
00487
00488 opac.OpacStack[i-oxy.i2d+opac.iopo2d] = crs;
00489 }
00490 opac.nOpacTot += iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon - oxy.i2d + 1;
00491
00492
00493
00494
00495 opac.ipOpMgEx = opac.nOpacTot + 1;
00496 for( i=opac.ipmgex-1; i < iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon; i++ )
00497 {
00498 opac.OpacStack[i-opac.ipmgex+opac.ipOpMgEx] =
00499 (0.2602325880970085 +
00500 445.8558249365131*exp(-rfield.AnuOrg[i]/0.1009243952792674))*
00501 1e-18;
00502 }
00503 opac.nOpacTot += iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon - opac.ipmgex + 1;
00504
00505 ASSERT( opac.nOpacTot < ndimOpacityStack );
00506
00507
00508
00509 OpacityCreatePowerLaw(opac.ica2ex[0],opac.ica2ex[1],4e-18,1.,&opac.ica2op);
00510
00511 if( trace.lgTrace )
00512 {
00513 fprintf( ioQQQ,
00514 " OpacityCreateAll return OK, number of opacity cells used in OPSC= %ld and OPSV is dimensioned %ld\n",
00515 opac.nOpacTot, ndimOpacityStack );
00516 }
00517
00518
00519
00520 if( opac.lgCompileOpac )
00521 {
00522 fprintf( ioQQQ, "The COMPILE OPACITIES command is currently not supported\n" );
00523 cdEXIT(EXIT_FAILURE);
00524 }
00525
00526 if( lgRealloc )
00527 fprintf(ioQQQ,
00528 " Please consider revising ndimOpacityStack in OpacityCreateAll, a total of %li cells were needed.\n\n" , opac.nOpacTot);
00529 return;
00530 }
00531
00532 STATIC void OpacityCreatePowerLaw(
00533
00534 long int ilo,
00535
00536 long int ihi,
00537
00538 double cross,
00539
00540 double s,
00541
00542 long int *ip)
00543 {
00544 long int i;
00545 double thres;
00546
00547 DEBUG_ENTRY( "OpacityCreatePowerLaw()" );
00548
00549
00550 ASSERT( cross > 0. );
00551
00552
00553 *ip = opac.nOpacTot + 1;
00554 ASSERT( *ip > 0 );
00555 ASSERT( ilo > 0 );
00556 thres = rfield.anu[ilo-1];
00557
00558 if( opac.nOpacTot + ihi - ilo + 1 > ndimOpacityStack )
00559 opacity_more_memory();
00560
00561 for( i=ilo-1; i < ihi; i++ )
00562 {
00563 opac.OpacStack[i-ilo+*ip] = cross*pow(rfield.anu[i]/thres,-s);
00564 }
00565
00566 opac.nOpacTot += ihi - ilo + 1;
00567 return;
00568 }
00569
00570
00571 STATIC void OpacityCreateReilMan(long int low,
00572 long int ihi,
00573 const realnum cross[],
00574 long int ncross,
00575 long int *ipop,
00576 const char *chLabl)
00577 {
00578 long int i,
00579 ics,
00580 j,
00581 ncr;
00582
00583 const int NOP = 100;
00584 realnum cs[NOP],
00585 en[NOP],
00586 slope;
00587
00588 DEBUG_ENTRY( "OpacityCreateReilMan()" );
00589
00590
00591
00592
00593
00594 *ipop = opac.nOpacTot + 1;
00595 ASSERT( *ipop > 0 );
00596
00597 ncr = ncross/2;
00598 if( ncr > NOP )
00599 {
00600 fprintf( ioQQQ, " Too many opacities were entered into OpacityCreateReilMan. Increase the value of NOP.\n" );
00601 fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl );
00602 cdEXIT(EXIT_FAILURE);
00603 }
00604
00605
00606
00607
00608 for( i=0; i < ncr; i++ )
00609 {
00610 en[i] = cross[i*2]/13.6f;
00611 cs[i] = cross[(i+1)*2-1]*1e-18f;
00612 }
00613
00614 ASSERT( low>0 );
00615 if( en[0] > rfield.anu[low-1] )
00616 {
00617 fprintf( ioQQQ,
00618 " OpacityCreateReilMan: The entered opacity energy bandwidth is not large enough (low fail).\n" );
00619 fprintf( ioQQQ,
00620 " The desired energy (Ryd) was%12.5eeV and the lowest entered in the array was%12.5e eV\n",
00621 rfield.anu[low-1]*EVRYD, en[0]*EVRYD );
00622
00623 fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl );
00624 fprintf( ioQQQ, " The original energy (eV) and cross section (mb) arrays follow:\n" );
00625 fprintf( ioQQQ, " " );
00626
00627 for( i=0; i < ncross; i++ )
00628 {
00629 fprintf( ioQQQ, "%11.4e", cross[i] );
00630 }
00631
00632 fprintf( ioQQQ, "\n" );
00633 cdEXIT(EXIT_FAILURE);
00634 }
00635
00636 slope = (cs[1] - cs[0])/(en[1] - en[0]);
00637 ics = 1;
00638
00639 if( opac.nOpacTot + ihi - low + 1 > ndimOpacityStack )
00640 opacity_more_memory();
00641
00642
00643 for( i=low-1; i < ihi; i++ )
00644 {
00645 if( rfield.anu[i] > en[ics-1] && rfield.anu[i] <= en[ics] )
00646 {
00647 opac.OpacStack[i-low+*ipop] = cs[ics-1] + slope*(rfield.anu[i] -
00648 en[ics-1]);
00649 }
00650
00651 else
00652 {
00653 ics += 1;
00654 if( ics + 1 > ncr )
00655 {
00656 fprintf( ioQQQ, " OpacityCreateReilMan: The entered opacity energy bandwidth is not large enough (high fail).\n" );
00657 fprintf( ioQQQ, " The entered energy was %10.2eeV and the highest in the array was %10.2eeV\n",
00658 rfield.anu[i]*13.6, en[ncr-1]*13.6 );
00659 fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl
00660 );
00661 fprintf( ioQQQ, " The lowest energy enterd in the array was%10.2e eV\n",
00662 en[0]*13.65 );
00663 fprintf( ioQQQ, " The highest energy ever needed would be%10.2eeV\n",
00664 rfield.anu[ihi-1]*13.6 );
00665 fprintf( ioQQQ, " The lowest energy needed was%10.2eeV\n",
00666 rfield.anu[low-1]*13.6 );
00667 cdEXIT(EXIT_FAILURE);
00668 }
00669
00670 slope = (cs[ics] - cs[ics-1])/(en[ics] - en[ics-1]);
00671 if( rfield.anu[i] > en[ics-1] && rfield.anu[i] <= en[ics] )
00672 {
00673 opac.OpacStack[i-low+*ipop] = cs[ics-1] + slope*(rfield.anu[i] -
00674 en[ics-1]);
00675 }
00676 else
00677 {
00678 ASSERT( i > 0);
00679 fprintf( ioQQQ, " Internal logical error in OpacityCreateReilMan.\n" );
00680 fprintf( ioQQQ, " The desired energy (%10.2eeV), I=%5ld, is not within the next energy bound%10.2e%10.2e\n",
00681 rfield.anu[i]*13.6, i, en[ics-1], en[ics] );
00682
00683 fprintf( ioQQQ, " The previous energy (eV) was%10.2e\n",
00684 rfield.anu[i-1]*13.6 );
00685
00686 fprintf( ioQQQ, " Here comes the energy array. ICS=%4ld\n",
00687 ics );
00688
00689 for( j=0; j < ncr; j++ )
00690 {
00691 fprintf( ioQQQ, "%10.2e", en[j] );
00692 }
00693 fprintf( ioQQQ, "\n" );
00694
00695 fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl );
00696 cdEXIT(EXIT_FAILURE);
00697 }
00698 }
00699 }
00700
00701
00702 opac.nOpacTot += ihi - low + 1;
00703 return;
00704 }
00705
00706
00707
00708 STATIC double ofit(double e,
00709 realnum opart[])
00710 {
00711 double otot,
00712 q,
00713 x;
00714
00715 static const double y[7][5] = {
00716 {8.915,3995.,3.242,10.44,0.0},
00717 {11.31,1498.,5.27,7.319,0.0},
00718 {10.5,1.059e05,1.263,13.04,0.0},
00719 {19.49,48.47,8.806,5.983,0.0},
00720 {50.,4.244e04,0.1913,7.012,4.454e-02},
00721 {110.5,0.1588,148.3,-3.38,3.589e-02},
00722 {177.4,32.37,381.2,1.083,0.0}
00723 };
00724 static const double eth[7]={13.62,16.94,18.79,28.48,50.,110.5,538.};
00725 static const long l[7]={1,1,1,0,1,1,0};
00726
00727 DEBUG_ENTRY( "ofit()" );
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740 otot = 0.0;
00741
00742 for( int i=0; i < 7; i++ )
00743 {
00744 opart[i] = 0.0;
00745 }
00746
00747 for( int i=0; i < 7; i++ )
00748 {
00749 if( e >= eth[i] )
00750 {
00751
00752 ASSERT( i < 7 );
00753 q = 5.5 - 0.5*y[i][3] + l[i];
00754
00755 x = e/y[i][0];
00756
00757 opart[i] = (realnum)(y[i][1]*(POW2(x - 1.0) + POW2(y[i][4]))/
00758 pow(x,q)/pow(1.0 + sqrt(x/y[i][2]),y[i][3]));
00759
00760 otot += opart[i];
00761 }
00762 }
00763 return(otot);
00764 }
00765
00766
00767
00768
00769 STATIC void OpacityCreate1Element(
00770
00771 long int nelem)
00772 {
00773 long int ihi,
00774 ip,
00775 ipop,
00776 limit,
00777 low,
00778 need,
00779 nelec,
00780 ion,
00781 nshell;
00782 double cs;
00783 double energy;
00784
00785 DEBUG_ENTRY( "OpacityCreate1Element()" );
00786
00787
00788 ASSERT( nelem >= 2 );
00789 ASSERT( nelem < LIMELM );
00790
00791
00792
00793 for( ion=0; ion < nelem; ion++ )
00794 {
00795
00796
00797 for( ip=0; ip < rfield.nupper; ip++ )
00798 {
00799 opac.opacity_abs[ip] = 0.;
00800 }
00801
00802
00803 nelec = nelem+1 - ion;
00804
00805
00806 for( nshell=0; nshell < Heavy.nsShells[nelem][ion]; nshell++ )
00807 {
00808
00809 opac.ipElement[nelem][ion][nshell][2] = opac.nOpacTot + 1;
00810
00811
00812 limit = opac.ipElement[nelem][ion][nshell][1];
00813
00814
00815 need = limit - opac.ipElement[nelem][ion][nshell][0] + 1;
00816
00817
00818 if( opac.nOpacTot + need > ndimOpacityStack )
00819 opacity_more_memory();
00820
00821
00822 low = opac.ipElement[nelem][ion][nshell][0];
00823 ihi = opac.ipElement[nelem][ion][nshell][1];
00824 ipop = opac.ipElement[nelem][ion][nshell][2];
00825
00826
00827
00828 ASSERT( low <= ihi || low<5 );
00829
00830
00831 for( ip=low-1; ip < ihi; ip++ )
00832 {
00833
00834 energy = MAX2(rfield.AnuOrg[ip]*EVRYD ,
00835 t_ADfA::Inst().ph1(nshell,nelec-1,nelem,0));
00836
00837
00838 cs = t_ADfA::Inst().phfit(nelem+1,nelec,nshell+1,energy);
00839
00840
00841
00842 opac.OpacStack[ip-low+ipop] = cs*1e-18;
00843
00844
00845 opac.opacity_abs[ip] += cs;
00846 }
00847
00848 opac.nOpacTot += ihi - low + 1;
00849
00850
00851 if( save.lgPunPoint )
00852 {
00853 fprintf( save.ipPoint, "%3ld%3ld%3ld%10.2e%10.2e%10.2e%10.2e\n",
00854 nelem, ion, nshell, rfield.anu[low-1], rfield.anu[ihi-1],
00855 opac.OpacStack[ipop-1], opac.OpacStack[ihi-low+ipop-1] );
00856 }
00857 }
00858
00859 ASSERT( Heavy.nsShells[nelem][ion] >= 1 );
00860
00861 for( ip=opac.ipElement[nelem][ion][Heavy.nsShells[nelem][ion]-1][0]-1;
00862 ip < continuum.KshellLimit; ip++ )
00863 {
00864 ASSERT( opac.opacity_abs[ip] > 0. );
00865 }
00866
00867 }
00868 return;
00869 }
00870
00871
00872 STATIC void opacity_more_memory(void)
00873 {
00874
00875 DEBUG_ENTRY( "opacity_more_memory()" );
00876
00877
00878 ndimOpacityStack *= 2;
00879 opac.OpacStack = (double *)REALLOC( opac.OpacStack , (size_t)ndimOpacityStack*sizeof(double) );
00880 fprintf( ioQQQ, " NOTE OpacityCreate1Element needed more opacity cells than ndimOpacityStack, please consider increasing it.\n" );
00881 fprintf( ioQQQ, " NOTE OpacityCreate1Element doubled memory allocation to %li.\n",ndimOpacityStack );
00882 lgRealloc = true;
00883 return;
00884 }
00885
00886
00887 STATIC double Opacity_iso_photo_cs(
00888
00889 double EgammaRyd ,
00890
00891 long ipISO ,
00892
00893 long nelem ,
00894
00895 long index )
00896 {
00897 double crs=-DBL_MAX;
00898
00899 DEBUG_ENTRY( "Opacity_iso_photo_cs()" );
00900
00901 if( ipISO==ipH_LIKE )
00902 {
00903 if( index==0 )
00904 {
00905
00906
00907 double EgammaEV = MAX2(EgammaRyd*(realnum)EVRYD , t_ADfA::Inst().ph1(0,0,nelem,0));
00908 crs = t_ADfA::Inst().phfit(nelem+1,1,1,EgammaEV)* 1e-18;
00909
00910 ASSERT( crs > 0. && crs < 1e-10 );
00911 }
00912 else if( index < iso_sp[ipISO][nelem].numLevels_max - iso_sp[ipISO][nelem].nCollapsed_max )
00913 {
00914
00915 double photon = MAX2( EgammaRyd/iso_sp[ipISO][nelem].fb[index].xIsoLevNIonRyd, 1. + FLT_EPSILON*2. );
00916
00917 crs = H_photo_cs( photon , N_(index), L_(index), nelem+1 );
00918
00919 ASSERT( crs > 0. && crs < 1e-10 );
00920 }
00921 else if( N_(index) <= NHYDRO_MAX_LEVEL )
00922 {
00923
00924
00925
00926
00927
00928
00929 EgammaRyd = MAX2( EgammaRyd , iso_sp[ipISO][nelem].fb[index].xIsoLevNIonRyd*1.001f );
00930
00931 crs = t_ADfA::Inst().hpfit(nelem+1,N_(index),EgammaRyd*EVRYD);
00932
00933 ASSERT( crs > 0. && crs < 1e-10 );
00934 }
00935 else
00936 {
00937
00938 double photon = MAX2( EgammaRyd/iso_sp[ipISO][nelem].fb[index].xIsoLevNIonRyd, 1. + FLT_EPSILON*2. );
00939
00940
00941
00942
00943 crs = H_photo_cs( photon , N_(index), N_(index)-1, nelem+1 );
00944
00945
00946 ASSERT( crs > 0. && crs < 1e-10 );
00947 }
00948 }
00949 else if( ipISO==ipHE_LIKE )
00950 {
00951 EgammaRyd = MAX2( EgammaRyd , iso_sp[ipISO][nelem].fb[index].xIsoLevNIonRyd);
00952
00953 if( index >= iso_sp[ipHE_LIKE][nelem].numLevels_max - iso_sp[ipHE_LIKE][nelem].nCollapsed_max )
00954 {
00955 long int nup = iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max + index + 1 -
00956 (iso_sp[ipHE_LIKE][nelem].numLevels_max - iso_sp[ipHE_LIKE][nelem].nCollapsed_max);
00957
00958
00959
00960 crs = t_ADfA::Inst().hpfit(nelem,nup ,EgammaRyd*EVRYD);
00961
00962 ASSERT(
00963 (EgammaRyd < iso_sp[ipISO][nelem].fb[index].xIsoLevNIonRyd*1.02) ||
00964 (crs > 0. && crs < 1e-10) );
00965 }
00966 else
00967 {
00968 long n = N_(index);
00969 long l = L_(index);
00970 long S = S_(index);
00971
00972
00973
00974
00975
00976 crs = He_cross_section( EgammaRyd, iso_sp[ipISO][nelem].fb[index].xIsoLevNIonRyd, n, l, S, nelem );
00977
00978
00979 ASSERT( crs > 0. && crs < 1e-10 );
00980 }
00981 }
00982 else
00983 TotalInsanity();
00984 return(crs);
00985 }
00986
00987
00988 static const int NCRS = 33;
00989
00990 STATIC double hmiopc(double freq)
00991 {
00992 double energy,
00993 hmiopc_v,
00994 x,
00995 y;
00996 static double y2[NCRS];
00997 static double crs[NCRS]={0.,0.124,0.398,0.708,1.054,1.437,1.805,
00998 2.176,2.518,2.842,3.126,3.377,3.580,3.741,3.851,3.913,3.925,
00999 3.887,3.805,3.676,3.511,3.306,3.071,2.810,2.523,2.219,1.898,
01000 1.567,1.233,.912,.629,.39,.19};
01001 static double ener[NCRS]={0.,0.001459,0.003296,0.005256,0.007351,
01002 0.009595,0.01201,0.01460,0.01741,0.02044,0.02375,0.02735,0.03129,
01003 0.03563,0.04043,0.04576,0.05171,0.05841,0.06601,0.07469,0.08470,
01004 0.09638,0.1102,0.1268,0.1470,0.1723,0.2049,0.2483,0.3090,0.4001,
01005 0.5520,0.8557,1.7669};
01006 static bool lgFirst = true;
01007
01008 DEBUG_ENTRY( "hmiopc()" );
01009
01010
01011
01012
01013
01014
01015
01016 if( lgFirst )
01017 {
01018
01019 spline(ener,crs,NCRS,2e31,2e31,y2);
01020 lgFirst = false;
01021 }
01022
01023 energy = freq - 0.05552;
01024 if( energy < ener[0] || energy > ener[NCRS-1] )
01025 {
01026 hmiopc_v = 0.;
01027 }
01028 else
01029 {
01030 x = energy;
01031 splint(ener,crs,y2,NCRS,x,&y);
01032 hmiopc_v = y*1e-17;
01033 }
01034 return( hmiopc_v );
01035 }
01036
01037
01038 STATIC double rayleh(double ener)
01039 {
01040 double rayleh_v;
01041
01042 DEBUG_ENTRY( "rayleh()" );
01043
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055 if( ener < 0.05 )
01056 {
01057 rayleh_v = (8.41e-25*powi(ener,4) + 3.37e-24*powi(ener,6))*
01058 hydro.DampOnFac;
01059 }
01060
01061 else if( ener < 0.646 )
01062 {
01063 rayleh_v = (8.41e-25*powi(ener,4) + 3.37e-24*powi(ener,6) +
01064 4.71e-22*powi(ener,14))*hydro.DampOnFac;
01065 }
01066
01067 else if( ener >= 0.646 && ener < 1.0 )
01068 {
01069 rayleh_v = fabs(0.74959-ener);
01070 rayleh_v = 1.788e5/POW2(FR1RYD*MAX2(0.001,rayleh_v));
01071
01072 rayleh_v = MAX2(rayleh_v,1e-24)*hydro.DampOnFac;
01073 }
01074
01075 else
01076 {
01077 rayleh_v = 0.;
01078 }
01079 return( rayleh_v );
01080 }