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