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 "path.h"
00028 #include "atmdat.h"
00029 #include "punch.h"
00030 #include "grains.h"
00031 #include "thirdparty.h"
00032 #include "hydro_bauman.h"
00033 #include "opacity.h"
00034 #include "helike_recom.h"
00035
00036 static const int NCSH2P = 10;
00037
00038
00039 static const int NMAGIC = 7;
00040
00041
00042 static const char MAGIC[NMAGIC] = "010118";
00043
00044
00045 static long int ndimOpacityStack = 1700000L;
00046
00047
00048
00049
00050
00051
00052 static double Yan_H2_CS( double energy_ryd )
00053 {
00054
00055 double energy_keV;
00056 double cross_section;
00057 double x;
00058 double xsqrt , x15 , x2;
00059 double energy = energy_ryd * EVRYD;
00060
00061 DEBUG_ENTRY( "Yan_H2_CS()" );
00062
00063
00064 x = energy / 15.4;
00065 energy_keV = energy/1000.0;
00066 xsqrt = sqrt(x);
00067 x15 = x * xsqrt;
00068 x2 = x*x;
00069
00070 if( energy < 15.4 )
00071 {
00072 cross_section = 0.;
00073 }
00074
00075 else if(energy >= 15.4 && energy < 18.)
00076 {
00077 cross_section = 1e7 * (1 - 197.448/xsqrt + 438.823/x - 260.481/x15 + 17.915/x2);
00078
00079 cross_section = MAX2( 0. , cross_section );
00080 }
00081
00082 else if(energy >= 18. && energy <= 30.)
00083 {
00084 cross_section = (-145.528 +351.394*xsqrt - 274.294*x + 74.320*x15)/pow(energy_keV,3.5);
00085 }
00086
00087 else if(energy > 30. && energy <= 85.)
00088 {
00089 cross_section = (65.304 - 91.762*xsqrt + 51.778*x - 9.364*x15)/pow(energy_keV,3.5);
00090 }
00091
00092
00093 else
00094 {
00095 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);
00096 }
00097
00098 DEBUG_EXIT( "Yan_H2_CS()" );
00099
00100 return( cross_section * 1e-24 );
00101 }
00102
00103
00104 static void OpacityCreate1Element(long int nelem);
00105
00106
00107 static void opacity_more_memory(void);
00108
00109
00110 static double hmiopc(double freq);
00111
00112
00113 static double rayleh(double ener);
00114
00115
00116 static double Opacity_iso_photo_cs( float energy , long ipISO , long nelem , long n );
00117
00118
00119 static void OpacityCreateReilMan(long int low,
00120 long int ihi,
00121 float cross[],
00122 long int ncross,
00123 long int *ipop,
00124 const char *chLabl);
00125
00126 static bool lgRealloc = false;
00127
00128
00129 static void OpacityCreatePowerLaw(
00130
00131 long int ilo,
00132
00133 long int ihi,
00134
00135 double cross,
00136
00137 double s,
00138
00139 long int *ip);
00140
00141
00142 static double ofit(double e,
00143 float opart[]);
00144
00145
00146 static void OpacityValenceRescale(
00147
00148 long int nelem ,
00149
00150 double scale )
00151 {
00152
00153 long int ion , nshell , low , ihi , ipop , ip;
00154
00155 DEBUG_ENTRY( "OpacityValenceRescale()" );
00156
00157
00158
00159
00160 if( !dense.lgElmtOn[nelem] )
00161 {
00162 DEBUG_EXIT( "OpacityValenceRescale()" );
00163 return;
00164 }
00165
00166
00167 ASSERT( scale >= 0. );
00168
00169 ion = 0;
00170
00171 nshell = Heavy.nsShells[nelem][ion] - 1;
00172
00173
00174 low = opac.ipElement[nelem][ion][nshell][0];
00175 ihi = opac.ipElement[nelem][ion][nshell][1];
00176 ipop = opac.ipElement[nelem][ion][nshell][2];
00177
00178
00179 for( ip=low-1; ip < ihi; ip++ )
00180 {
00181 opac.OpacStack[ip-low+ipop] *= scale;
00182 }
00183
00184 DEBUG_EXIT( "OpacityValenceRescale()" );
00185 return;
00186 }
00187
00188 void OpacityCreateAll(void)
00189 {
00190 long int i,
00191 ipISO ,
00192 n ,
00193 need ,
00194 nelem;
00195
00196 float opart[7];
00197
00198 double crs,
00199 dx,
00200 eps,
00201 thom,
00202 thres,
00203 x;
00204
00205
00206 char chFileName[FILENAME_PATH_LENGTH_2];
00207
00208 FILE *ioOPAC;
00209
00210
00211
00212
00213
00214
00215 static float csh2p[NCSH2P]={6.75f,0.24f,8.68f,2.5f,10.54f,7.1f,12.46f,
00216 6.0f,14.28f,2.7f};
00217
00218 DEBUG_ENTRY( "OpacityCreateAll()" );
00219
00220
00221
00222
00223
00224
00225 GrainsInit();
00226
00227
00228
00229
00230 if( lgOpacMalloced )
00231 {
00232
00233 if( trace.lgTrace )
00234 {
00235 fprintf( ioQQQ, " OpacityCreateAll called but NOT evaluated since already done.\n" );
00236 }
00237
00238 DEBUG_EXIT( "OpacityCreateAll()" );
00239 return;
00240 }
00241
00242 lgOpacMalloced = true;
00243
00244
00245 opac.OpacStack = (double*)MALLOC((size_t)ndimOpacityStack*sizeof(double));
00246
00247
00248 if( opac.lgUseFileOpac )
00249 {
00250
00251
00252 strcpy( chFileName , "opacity.opc" );
00253 ioOPAC = fopen( chFileName , "rb" );
00254 if( ioOPAC == NULL )
00255 {
00256
00257 if( lgDataPathSet )
00258 {
00259
00260 strcpy( chFileName , chDataPath );
00261 strcat( chFileName , "opacity.opc" );
00262
00263 ioOPAC = fopen( chFileName , "rb" );
00264 if( ioOPAC == NULL )
00265 {
00266
00267 opac.lgOpacExist = false;
00268 }
00269 else
00270 {
00271
00272 opac.lgOpacExist = true;
00273 }
00274 }
00275 else
00276 {
00277 opac.lgOpacExist = false;
00278 }
00279 }
00280 else
00281 {
00282 opac.lgOpacExist = true;
00283 }
00284
00285
00286 if( !opac.lgCompileOpac && opac.lgOpacExist && (ioOPAC!=NULL) )
00287 {
00288 char chMagic[NMAGIC];
00289
00290
00291
00292 n = (long)fread( chMagic , NMAGIC , sizeof(char), ioOPAC );
00293 if( strcmp( chMagic , MAGIC ) != 0 )
00294 {
00295 fprintf( ioQQQ," The starting magic number in the compiled opacity file is incorrect.\n");
00296 fprintf( ioQQQ," Please recompile opacities, or delete opacity.opc in the data directory.\n");
00297 fprintf( ioQQQ," The magic number was %s\n", MAGIC );
00298
00299 path_not_set();
00300 puts( "[Stop in OpacityCreateAll]" );
00301 cdEXIT(EXIT_FAILURE);
00302 }
00303
00304
00305 n = (long)fread( opac.OpacStack , 1, sizeof(opac.OpacStack), ioOPAC );
00306 if( ((unsigned)n-sizeof(opac.OpacStack)) != 0 )
00307 {
00308 fprintf( ioQQQ, " problem trying to read opacity.opc\n" );
00309 fprintf( ioQQQ, " I expected to read %li words, but fread returned only %li\n",
00310 (long)sizeof(opac.OpacStack),n);
00311 fprintf( ioQQQ, " Try recompiling the opacities with the COMPILE OPACITY command,\n" );
00312 fprintf( ioQQQ, " or tell the code not to use the compiled opacities with the NO FILE OPACITY command.\n" );
00313
00314
00315 if( feof(ioOPAC ) )
00316 {
00317 fprintf( ioQQQ, " end of file hit\n" );
00318 }
00319
00320 if( ferror(ioOPAC ) )
00321 {
00322 fprintf( ioQQQ, " end of file hit\n" );
00323 }
00324 puts( "[Stop in OpacityCreateAll]" );
00325 cdEXIT(EXIT_FAILURE);
00326 }
00327
00328
00329
00330 n = 0;
00331
00332 n -= (long)(fread( iso.ipOpac[ipH_LIKE], 1,sizeof(iso.ipOpac[ipH_LIKE]),
00333 ioOPAC ) - sizeof(iso.ipOpac[ipH_LIKE]));
00334
00335 n -= (long)(fread( &opac.ipRayScat,1,sizeof(long) ,
00336 ioOPAC ) - sizeof(long));
00337
00338 n -= (long)(fread( &opac.iopcom, 1,sizeof(long) ,
00339 ioOPAC ) - sizeof(long));
00340
00341 n -= (long)(fread( &opac.ippr, 1,sizeof(long) ,
00342 ioOPAC ) - sizeof(long));
00343
00344 n -= (long)(fread( &opac.ioppr, 1,sizeof(long) ,
00345 ioOPAC ) - sizeof(long));
00346
00347 n -= (long)(fread( &opac.ipBrems, 1,sizeof(long) ,
00348 ioOPAC ) - sizeof(long));
00349
00350 n -= (long)(fread( &opac.iphmra, 1,sizeof(long) ,
00351 ioOPAC ) - sizeof(long));
00352
00353 n -= (long)(fread( &opac.iphmop, 1,sizeof(long) ,
00354 ioOPAC ) - sizeof(long));
00355
00356 n -= (long)(fread( opac.ih2pnt, 1,sizeof(opac.ih2pnt),
00357 ioOPAC ) - sizeof(opac.ih2pnt));
00358
00359 n -= (long)(fread( &opac.ih2pof, 1,sizeof(long) ,
00360 ioOPAC ) - sizeof(long));
00361
00362 n -= (long)(fread( opac.iophe1, 1,sizeof(opac.iophe1),
00363 ioOPAC ) - sizeof(opac.iophe1));
00364
00365 n -= (long)(fread( &opac.ioptri, 1,sizeof(long) ,
00366 ioOPAC ) - sizeof(long));
00367
00368 n -= (long)(fread( opac.ipElement, 1,sizeof(opac.ipElement),
00369 ioOPAC ) - sizeof(opac.ipElement));
00370
00371 n -= (long)(fread( opac.in1, 1,sizeof(opac.in1),
00372 ioOPAC ) - sizeof(opac.in1));
00373
00374 n -= (long)(fread( opac.ipo3exc, 1,sizeof(opac.ipo3exc),
00375 ioOPAC ) - sizeof(opac.ipo3exc));
00376
00377 n -= (long)(fread( opac.ipo3exc3, 1,sizeof(opac.ipo3exc3),
00378 ioOPAC ) - sizeof(opac.ipo3exc3));
00379
00380 n -= (long)(fread( opac.ipo1exc, 1,sizeof(opac.ipo1exc),
00381 ioOPAC ) - sizeof(opac.ipo1exc));
00382
00383 n -= (long)(fread( &opac.iopo2d, 1,sizeof(long) ,
00384 ioOPAC ) - sizeof(long));
00385
00386 n -= (long)(fread( &opac.ipmgex, 1,sizeof(long) ,
00387 ioOPAC ) - sizeof(long));
00388
00389 n -= (long)(fread( &opac.ipOpMgEx, 1,sizeof(long) ,
00390 ioOPAC ) - sizeof(long));
00391
00392 n -= (long)(fread( opac.ica2ex, 1,sizeof(opac.ica2ex),
00393 ioOPAC ) - sizeof(opac.ica2ex));
00394
00395 n -= (long)(fread( &opac.ica2op, 1,sizeof(long) ,
00396 ioOPAC ) - sizeof(long));
00397
00398
00399 n -= (long)(fread( opac.tmn, 1,sizeof(opac.tmn),
00400 ioOPAC ) - sizeof(opac.tmn));
00401
00402
00403 if( n!= 0 )
00404 {
00405 fprintf( ioQQQ, " problem trying to read opacity.opc pointers\n" );
00406 fprintf( ioQQQ, " fread was short by %li\n",
00407 n);
00408 fprintf( ioQQQ, " Try recompiling the opacities with the COMPILE OPACITY command,\n" );
00409 fprintf( ioQQQ, " or tell the code not to use the compiled opacities with the NO FILE OPACITY command.\n" );
00410 puts( "[Stop in OpacityCreateAll]" );
00411 cdEXIT(EXIT_FAILURE);
00412 }
00413
00414
00415 if( opac.tmn[rfield.nupper-1] != rfield.AnuOrg[rfield.nupper-1] )
00416
00417 {
00418 fprintf( ioQQQ, " Energy grid in opacity file opacity.opc opacity.pnt do not agree with this version of the code - please recompile.\n" );
00419 fprintf( ioQQQ, " Values are:%12.4e%12.4e\n",
00420 rfield.AnuOrg[rfield.nupper-1], opac.tmn[rfield.nupper-1] );
00421 puts( "[Stop in OpacityCreateAll]" );
00422 cdEXIT(EXIT_FAILURE);
00423 }
00424
00425 if( trace.lgTrace )
00426 {
00427 fprintf( ioQQQ, " OpacityCreateAll called, values from file =%s=\n",chFileName );
00428 }
00429
00430
00431 opac.lgOpacExist = true;
00432
00433
00434 for( i=0; i < rfield.nupper; i++ )
00435 {
00436 opac.opacity_abs[i] = 0.;
00437 }
00438
00439 DEBUG_EXIT( "OpacityCreateAll()" );
00440 return;
00441 }
00442 else
00443 {
00444 opac.lgOpacExist = false;
00445 }
00446 }
00447 else
00448 {
00449
00450 opac.lgOpacExist = false;
00451 }
00452
00453 if( trace.lgTrace )
00454 {
00455 fprintf( ioQQQ, " OpacityCreateAll called, evaluating.\n" );
00456 }
00457
00458
00459 for( i=0; i < rfield.nupper; i++ )
00460 {
00461 opac.opacity_abs[i] = 0.;
00462 }
00463
00464
00465 opac.nOpacTot = 0;
00466
00467
00468 for( ipISO=ipH_LIKE; ipISO<=ipHE_LIKE; ++ipISO )
00469 {
00470 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00471 {
00472 if( dense.lgElmtOn[nelem] )
00473 {
00474 long int nupper;
00475
00476
00477
00478
00479 opac.ipElement[nelem][nelem-ipISO][0][2] = opac.nOpacTot + 1;
00480
00481
00482
00483 nupper = rfield.nupper;
00484 for( n=0; n < iso.numLevels_max[ipISO][nelem]; n++ )
00485 {
00486
00487 iso.ipOpac[ipISO][nelem][n] = opac.nOpacTot + 1;
00488
00489
00490
00491
00492 ASSERT( rfield.AnuOrg[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 0.94f *
00493 iso.xIsoLevNIonRyd[ipISO][nelem][n] );
00494
00495
00496 need = nupper - iso.ipIsoLevNIonCon[ipISO][nelem][n] + 1;
00497 ASSERT( need > 0 );
00498
00499 if( opac.nOpacTot + need > ndimOpacityStack )
00500 opacity_more_memory();
00501
00502 for( i=iso.ipIsoLevNIonCon[ipISO][nelem][n]-1; i < nupper; i++ )
00503 {
00504 opac.OpacStack[i-iso.ipIsoLevNIonCon[ipISO][nelem][n]+iso.ipOpac[ipISO][nelem][n]] =
00505 Opacity_iso_photo_cs( rfield.AnuOrg[i] , ipISO , nelem , n );
00506 }
00507
00508 opac.nOpacTot += need;
00509
00510 nupper = iso.ipIsoLevNIonCon[ipISO][nelem][0];
00511 }
00512 }
00513 }
00514 }
00515
00516
00517
00518 if( opac.nOpacTot + iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] + rfield.nupper > ndimOpacityStack )
00519 opacity_more_memory();
00520
00521
00522 opac.ipRayScat = opac.nOpacTot + 1;
00523 for( i=0; i < iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]; i++ )
00524 {
00525 opac.OpacStack[i-1+opac.ipRayScat] = rayleh(rfield.AnuOrg[i]);
00526 }
00527 opac.nOpacTot += iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s];
00528
00529
00530
00531
00532
00533
00534 thom = 6.65e-25;
00535 opac.iopcom = opac.nOpacTot + 1;
00536 for( i=0; i < opac.ipCKshell; i++ )
00537 {
00538 opac.OpacStack[i-1+opac.iopcom] = thom;
00539
00540
00541 }
00542
00543
00544
00545 for( i=opac.ipCKshell; i < rfield.nupper; i++ )
00546 {
00547 dx = rfield.AnuOrg[i]/3.7573e4;
00548
00549 opac.OpacStack[i-1+opac.iopcom] = thom*3.e0/4.e0*((1.e0 +
00550 dx)/POW3(dx)*(2.e0*dx*(1.e0 + dx)/(1.e0 + 2.e0*dx) - log(1.e0+
00551 2.e0*dx)) + 1.e0/2.e0/dx*log(1.e0+2.e0*dx) - (1.e0 + 3.e0*
00552 dx)/POW3(1.e0 + 2.e0*dx));
00553
00554
00555 }
00556 opac.nOpacTot += rfield.nupper - 1 + 1;
00557
00558
00559
00560
00561
00562 if( opac.nOpacTot + 3*rfield.nupper - opac.ippr + iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] - hmi.iphmin + 2 > ndimOpacityStack )
00563 opacity_more_memory();
00564
00565
00566 opac.ioppr = opac.nOpacTot + 1;
00567 for( i=opac.ippr-1; i < rfield.nupper; i++ )
00568 {
00569
00570
00571
00572
00573 x = rfield.AnuOrg[i]/7.512e4*2.;
00574
00575 opac.OpacStack[i-opac.ippr+opac.ioppr] = 5.793e-28*
00576 POW2((-0.46737118 + x*(0.349255416 + x*0.002179893))/(1. +
00577 x*(0.130471301 + x*0.000524906)));
00578
00579
00580 }
00581 opac.nOpacTot += rfield.nupper - opac.ippr + 1;
00582
00583
00584 opac.ipBrems = opac.nOpacTot + 1;
00585
00586 for( i=0; i < rfield.nupper; i++ )
00587 {
00588
00589
00590 opac.OpacStack[i-1+opac.ipBrems] =
00591
00592
00593
00594 1.03680e-8/POW3(rfield.AnuOrg[i]);
00595 }
00596 opac.nOpacTot += rfield.nupper - 1 + 1;
00597
00598 opac.iphmra = opac.nOpacTot + 1;
00599 for( i=0; i < rfield.nupper; i++ )
00600 {
00601
00602 opac.OpacStack[i-1+opac.iphmra] = 0.1175*rfield.anusqr[i];
00603 }
00604 opac.nOpacTot += rfield.nupper - 1 + 1;
00605
00606 opac.iphmop = opac.nOpacTot + 1;
00607 for( i=hmi.iphmin-1; i < iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]; i++ )
00608 {
00609
00610 opac.OpacStack[i-hmi.iphmin+opac.iphmop] =
00611 hmiopc(rfield.AnuOrg[i]);
00612 }
00613 opac.nOpacTot += iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] - hmi.iphmin + 1;
00614
00615 opac.ipH2_photo_opac_offset = opac.nOpacTot + 1;
00616 for( i=opac.ipH2_photo_thresh-1; i < rfield.nupper; i++ )
00617 {
00618
00619 opac.OpacStack[i-opac.ipH2_photo_thresh + opac.ipH2_photo_opac_offset] =
00620 Yan_H2_CS(rfield.AnuOrg[i]);
00621
00622
00623
00624
00625 }
00626 opac.nOpacTot += rfield.nupper - opac.ipH2_photo_thresh + 1;
00627
00628
00629
00630
00631 OpacityCreateReilMan(opac.ih2pnt[0],opac.ih2pnt[1],csh2p,NCSH2P,&opac.ih2pof,
00632 "H2+ ");
00633
00634
00635
00636 if( opac.nOpacTot + rfield.nupper - iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] + 1 > ndimOpacityStack )
00637 opacity_more_memory();
00638
00639
00640 opac.iophe1[0] = opac.nOpacTot + 1;
00641 opac.ipElement[ipHELIUM][0][0][2] = opac.iophe1[0];
00642 for( i=iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]-1; i < rfield.nupper; i++ )
00643 {
00644 crs = t_ADfA::Inst().phfit(2,2,1,rfield.AnuOrg[i]*EVRYD);
00645 opac.OpacStack[i-iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]+opac.iophe1[0]] =
00646 crs*1e-18;
00647 }
00648 opac.nOpacTot += rfield.nupper - iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] + 1;
00649
00650
00651
00652
00653
00654
00655 for( nelem=2; nelem < LIMELM; nelem++ )
00656 {
00657 if( dense.lgElmtOn[nelem] )
00658 {
00659 OpacityCreate1Element(nelem);
00660 }
00661 }
00662
00663
00664
00665
00666
00667
00668 OpacityValenceRescale( ipPOTASSIUM , 5. );
00669
00670
00671
00672
00673
00674 OpacityCreatePowerLaw(opac.in1[0],opac.in1[1],9e-18,1.75,&opac.in1[2]);
00675
00676
00677
00678 if( dense.lgElmtOn[ipOXYGEN] && t_ADfA::Inst().get_version() == PHFIT96 )
00679 {
00680
00681
00682 if( opac.nOpacTot + opac.ipElement[ipOXYGEN][0][2][1] - opac.ipElement[ipOXYGEN][0][2][0] + 1 > ndimOpacityStack )
00683 opacity_more_memory();
00684
00685
00686 for( i=opac.ipElement[ipOXYGEN][0][2][0]-1; i < opac.ipElement[ipOXYGEN][0][2][1]; i++ )
00687 {
00688
00689 eps = rfield.AnuOrg[i]*EVRYD;
00690 crs = ofit(eps,opart);
00691
00692
00693 crs = opart[0];
00694 for( n=1; n < 6; n++ )
00695 {
00696
00697 crs += opart[n];
00698 }
00699
00700 crs *= 1e-18;
00701 opac.OpacStack[i-opac.ipElement[ipOXYGEN][0][2][0]+opac.ipElement[ipOXYGEN][0][2][2]] = crs;
00702 }
00703
00704
00705 opac.nOpacTot += opac.ipElement[ipOXYGEN][0][2][1] - opac.ipElement[ipOXYGEN][0][2][0] + 1;
00706 }
00707
00708
00709 OpacityCreatePowerLaw(opac.ipo1exc[0],opac.ipo1exc[1],4.64e-18,0.,&opac.ipo1exc[2]);
00710
00711
00712
00713 OpacityCreatePowerLaw(opac.ipo3exc[0],opac.ipo3exc[1],3.8e-18,0.,&opac.ipo3exc[2]);
00714
00715
00716 OpacityCreatePowerLaw(opac.ipo3exc3[0],opac.ipo3exc3[1],5.5e-18,0.01,
00717 &opac.ipo3exc3[2]);
00718
00719
00720
00721 if( opac.nOpacTot + iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s] - oxy.i2d + 1
00722 + iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - opac.ipmgex + 1 > ndimOpacityStack )
00723 opacity_more_memory();
00724
00725
00726 opac.iopo2d = opac.nOpacTot + 1;
00727 thres = rfield.AnuOrg[oxy.i2d-1];
00728 for( i=oxy.i2d-1; i < iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]; i++ )
00729 {
00730 crs = 3.85e-18*(4.4*pow(rfield.AnuOrg[i]/thres,-1.5) - 3.38*
00731 pow(rfield.AnuOrg[i]/thres,-2.5));
00732
00733 opac.OpacStack[i-oxy.i2d+opac.iopo2d] = crs;
00734 }
00735 opac.nOpacTot += iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s] - oxy.i2d + 1;
00736
00737
00738
00739
00740 opac.ipOpMgEx = opac.nOpacTot + 1;
00741 for( i=opac.ipmgex-1; i < iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]; i++ )
00742 {
00743 opac.OpacStack[i-opac.ipmgex+opac.ipOpMgEx] =
00744 (0.2602325880970085 +
00745 445.8558249365131*exp(-rfield.AnuOrg[i]/0.1009243952792674))*
00746 1e-18;
00747 }
00748 opac.nOpacTot += iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - opac.ipmgex + 1;
00749
00750 ASSERT( opac.nOpacTot < ndimOpacityStack );
00751
00752
00753
00754 OpacityCreatePowerLaw(opac.ica2ex[0],opac.ica2ex[1],4e-18,1.,&opac.ica2op);
00755
00756 if( trace.lgTrace )
00757 {
00758 fprintf( ioQQQ,
00759 " OpacityCreateAll return OK, number of opacity cells used in OPSC= %ld and OPSV is dimensioned %ld\n",
00760 opac.nOpacTot, ndimOpacityStack );
00761 }
00762
00763
00764
00765 if( opac.lgCompileOpac )
00766 {
00767 ioOPAC = fopen("opacity.opc","wb");
00768 if( ioOPAC == NULL )
00769 {
00770 fprintf( ioQQQ, " problem trying to open opacity.opc\n" );
00771 puts( "[Stop in OpacityCreateAll]" );
00772 cdEXIT(EXIT_FAILURE);
00773 }
00774
00775
00776 n = (long)fwrite( MAGIC , NMAGIC , sizeof(char), ioOPAC );
00777
00778
00779 n = (long)fwrite( opac.OpacStack , 1, sizeof(opac.OpacStack), ioOPAC );
00780 if( ((size_t)n-sizeof(opac.OpacStack)) != 0 )
00781 {
00782 fprintf( ioQQQ, " problem trying to write opacity.opc\n" );
00783 fprintf( ioQQQ, " I expected to write %li words, but fwrite returned only %li\n",
00784 (long)sizeof(opac.OpacStack),n);
00785 puts( "[Stop in OpacityCreateAll]" );
00786 cdEXIT(EXIT_FAILURE);
00787 }
00788
00789
00790
00791 n = 0;
00792
00793 n -= (long)( fwrite( iso.ipOpac[ipH_LIKE], 1,sizeof(iso.ipOpac[ipH_LIKE]),
00794 ioOPAC ) - sizeof(iso.ipOpac[ipH_LIKE]));
00795
00796 n -= (long)( fwrite( &opac.ipRayScat, 1,sizeof(long) ,
00797 ioOPAC ) - sizeof(long));
00798
00799 n -= (long)( fwrite( &opac.iopcom, 1,sizeof(long) ,
00800 ioOPAC ) - sizeof(long));
00801
00802 n -= (long)( fwrite( &opac.ippr, 1,sizeof(long) ,
00803 ioOPAC ) - sizeof(long));
00804
00805 n -= (long)( fwrite( &opac.ioppr, 1,sizeof(long) ,
00806 ioOPAC ) - sizeof(long));
00807
00808 n -= (long)( fwrite( &opac.ipBrems, 1,sizeof(long) ,
00809 ioOPAC ) - sizeof(long));
00810
00811 n -= (long)( fwrite( &opac.iphmra, 1,sizeof(long) ,
00812 ioOPAC ) - sizeof(long));
00813
00814 n -= (long)( fwrite( &opac.iphmop, 1,sizeof(long) ,
00815 ioOPAC ) - sizeof(long));
00816
00817 n -= (long)( fwrite( opac.ih2pnt, 1,sizeof(opac.ih2pnt),
00818 ioOPAC ) - sizeof(opac.ih2pnt));
00819
00820 n -= (long)( fwrite( &opac.ih2pof, 1,sizeof(long) ,
00821 ioOPAC ) - sizeof(long));
00822
00823 n -= (long)( fwrite( opac.iophe1, 1,sizeof(opac.iophe1),
00824 ioOPAC ) - sizeof(opac.iophe1));
00825
00826 n -= (long)( fwrite( &opac.ioptri, 1,sizeof(long) ,
00827 ioOPAC ) - sizeof(long));
00828
00829 n -= (long)( fwrite( opac.ipElement, 1,sizeof(opac.ipElement),
00830 ioOPAC ) - sizeof(opac.ipElement));
00831
00832 n -= (long)( fwrite( opac.in1, 1,sizeof(opac.in1),
00833 ioOPAC ) - sizeof(opac.in1));
00834
00835 n -= (long)( fwrite( opac.ipo3exc, 1,sizeof(opac.ipo3exc),
00836 ioOPAC ) - sizeof(opac.ipo3exc));
00837
00838 n -= (long)( fwrite( opac.ipo3exc3, 1,sizeof(opac.ipo3exc3),
00839 ioOPAC ) - sizeof(opac.ipo3exc3));
00840
00841 n -= (long)( fwrite( opac.ipo1exc, 1,sizeof(opac.ipo1exc),
00842 ioOPAC ) - sizeof(opac.ipo1exc));
00843
00844 n -= (long)( fwrite( &opac.iopo2d, 1,sizeof(long) ,
00845 ioOPAC ) - sizeof(long));
00846
00847 n -= (long)( fwrite( &opac.ipmgex, 1,sizeof(long) ,
00848 ioOPAC ) - sizeof(long));
00849
00850 n -= (long)( fwrite( &opac.ipOpMgEx, 1,sizeof(long) ,
00851 ioOPAC ) - sizeof(long));
00852
00853 n -= (long)( fwrite( opac.ica2ex, 1,sizeof(opac.ica2ex),
00854 ioOPAC ) - sizeof(opac.ica2ex));
00855
00856 n -= (long)( fwrite( &opac.ica2op, 1,sizeof(long) ,
00857 ioOPAC ) - sizeof(long));
00858
00859 n -= (long)( fwrite( rfield.AnuOrg, 1,sizeof(rfield.AnuOrg),
00860 ioOPAC ) - sizeof(rfield.AnuOrg));
00861
00862
00863 if( n!= 0 )
00864 {
00865 fprintf( ioQQQ, " problem trying to write opacity.opc pointers\n" );
00866 fprintf( ioQQQ, " fwrite was short by %li\n",
00867 n);
00868 puts( "[Stop in OpacityCreateAll]" );
00869 cdEXIT(EXIT_FAILURE);
00870 }
00871
00872 fclose(ioOPAC);
00873 fprintf( ioQQQ, "\n\nCompile opacity completed ok - stopping.\n" );
00874 fprintf( ioQQQ, "The file opacity.opc was created.\n" );
00875 fprintf( ioQQQ, "Make sure this file lies somewhere on the path.\n\n\n" );
00876 puts( "[Stop in OpacityCreateAll]" );
00877 cdEXIT(EXIT_FAILURE);
00878 }
00879
00880 if( lgRealloc )
00881 fprintf(ioQQQ,
00882 " Please consider revising ndimOpacityStack in OpacityCreateAll, a total of %li cells were needed.\n\n" , opac.nOpacTot);
00883
00884 DEBUG_EXIT( "OpacityCreateAll()" );
00885 return;
00886 }
00887
00888 static void OpacityCreatePowerLaw(
00889
00890 long int ilo,
00891
00892 long int ihi,
00893
00894 double cross,
00895
00896 double s,
00897
00898 long int *ip)
00899 {
00900 long int i;
00901 double thres;
00902
00903 DEBUG_ENTRY( "OpacityCreatePowerLaw()" );
00904
00905
00906 ASSERT( cross > 0. );
00907
00908
00909 *ip = opac.nOpacTot + 1;
00910 ASSERT( *ip > 0 );
00911 ASSERT( ilo > 0 );
00912 thres = rfield.anu[ilo-1];
00913
00914 if( opac.nOpacTot + ihi - ilo + 1 > ndimOpacityStack )
00915 opacity_more_memory();
00916
00917 for( i=ilo-1; i < ihi; i++ )
00918 {
00919 opac.OpacStack[i-ilo+*ip] = cross*pow(rfield.anu[i]/thres,-s);
00920 }
00921
00922 opac.nOpacTot += ihi - ilo + 1;
00923
00924 DEBUG_EXIT( "OpacityCreatePowerLaw()" );
00925 return;
00926 }
00927
00928
00929 static void OpacityCreateReilMan(long int low,
00930 long int ihi,
00931 float cross[],
00932 long int ncross,
00933 long int *ipop,
00934 const char *chLabl)
00935 {
00936 long int i,
00937 ics,
00938 j,
00939 ncr;
00940
00941 const int NOP = 100;
00942 float cs[NOP],
00943 en[NOP],
00944 slope;
00945
00946 DEBUG_ENTRY( "OpacityCreateReilMan()" );
00947
00948
00949
00950
00951
00952 *ipop = opac.nOpacTot + 1;
00953 ASSERT( *ipop > 0 );
00954
00955 ncr = ncross/2;
00956 if( ncr > NOP )
00957 {
00958 fprintf( ioQQQ, " Too many opacities were entered into OpacityCreateReilMan. Increase the value of NOP.\n" );
00959 fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl );
00960 puts( "[Stop in OpacityCreateReilMan]" );
00961 cdEXIT(EXIT_FAILURE);
00962 }
00963
00964
00965
00966
00967 for( i=0; i < ncr; i++ )
00968 {
00969 en[i] = (float)(cross[i*2]/13.6);
00970 cs[i] = (float)(cross[(i+1)*2-1]*1e-18);
00971 }
00972
00973 ASSERT( low>0 );
00974 if( en[0] > rfield.anu[low-1] )
00975 {
00976 fprintf( ioQQQ,
00977 " OpacityCreateReilMan: The entered opacity energy bandwidth is not large enough (low fail).\n" );
00978 fprintf( ioQQQ,
00979 " The desired energy (Ryd) was%12.5eeV and the lowest entered in the array was%12.5e eV\n",
00980 rfield.anu[low-1]*EVRYD, en[0]*EVRYD );
00981
00982 fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl );
00983 fprintf( ioQQQ, " The original energy (eV) and cross section (mb) arrays follow:\n" );
00984 fprintf( ioQQQ, " " );
00985
00986 for( i=0; i < ncross; i++ )
00987 {
00988 fprintf( ioQQQ, "%11.4e", cross[i] );
00989 }
00990
00991 fprintf( ioQQQ, "\n" );
00992 puts( "[Stop in OpacityCreateReilMan]" );
00993 cdEXIT(EXIT_FAILURE);
00994 }
00995
00996 slope = (cs[1] - cs[0])/(en[1] - en[0]);
00997 ics = 1;
00998
00999 if( opac.nOpacTot + ihi - low + 1 > ndimOpacityStack )
01000 opacity_more_memory();
01001
01002
01003 for( i=low-1; i < ihi; i++ )
01004 {
01005 if( rfield.anu[i] > en[ics-1] && rfield.anu[i] <= en[ics] )
01006 {
01007 opac.OpacStack[i-low+*ipop] = cs[ics-1] + slope*(rfield.anu[i] -
01008 en[ics-1]);
01009 }
01010
01011 else
01012 {
01013 ics += 1;
01014 if( ics + 1 > ncr )
01015 {
01016 fprintf( ioQQQ, " OpacityCreateReilMan: The entered opacity energy bandwidth is not large enough (high fail).\n" );
01017 fprintf( ioQQQ, " The entered energy was %10.2eeV and the highest in the array was %10.2eeV\n",
01018 rfield.anu[i]*13.6, en[ncr-1]*13.6 );
01019 fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl
01020 );
01021 fprintf( ioQQQ, " The lowest energy enterd in the array was%10.2e eV\n",
01022 en[0]*13.65 );
01023 fprintf( ioQQQ, " The highest energy ever needed would be%10.2eeV\n",
01024 rfield.anu[ihi-1]*13.6 );
01025 fprintf( ioQQQ, " The lowest energy needed was%10.2eeV\n",
01026 rfield.anu[low-1]*13.6 );
01027 puts( "[Stop in OpacityCreateReilMan]" );
01028 cdEXIT(EXIT_FAILURE);
01029 }
01030
01031 slope = (cs[ics] - cs[ics-1])/(en[ics] - en[ics-1]);
01032 if( rfield.anu[i] > en[ics-1] && rfield.anu[i] <= en[ics] )
01033 {
01034 opac.OpacStack[i-low+*ipop] = cs[ics-1] + slope*(rfield.anu[i] -
01035 en[ics-1]);
01036 }
01037 else
01038 {
01039 ASSERT( i > 0);
01040 fprintf( ioQQQ, " Internal logical error in OpacityCreateReilMan.\n" );
01041 fprintf( ioQQQ, " The desired energy (%10.2eeV), I=%5ld, is not within the next energy bound%10.2e%10.2e\n",
01042 rfield.anu[i]*13.6, i, en[ics-1], en[ics] );
01043
01044 fprintf( ioQQQ, " The previous energy (eV) was%10.2e\n",
01045 rfield.anu[i-1]*13.6 );
01046
01047 fprintf( ioQQQ, " Here comes the energy array. ICS=%4ld\n",
01048 ics );
01049
01050 for( j=0; j < ncr; j++ )
01051 {
01052 fprintf( ioQQQ, "%10.2e", en[j] );
01053 }
01054 fprintf( ioQQQ, "\n" );
01055
01056 fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl );
01057 puts( "[Stop in OpacityCreateReilMan]" );
01058 cdEXIT(EXIT_FAILURE);
01059 }
01060 }
01061 }
01062
01063
01064 opac.nOpacTot += ihi - low + 1;
01065
01066 DEBUG_EXIT( "OpacityCreateReilMan()" );
01067 return;
01068 }
01069
01070
01071
01072 static double ofit(double e,
01073 float opart[])
01074 {
01075 double otot,
01076 q,
01077 x;
01078
01079 static const double y[7][5] = {
01080 {8.915,3995.,3.242,10.44,0.0},
01081 {11.31,1498.,5.27,7.319,0.0},
01082 {10.5,1.059e05,1.263,13.04,0.0},
01083 {19.49,48.47,8.806,5.983,0.0},
01084 {50.,4.244e04,0.1913,7.012,4.454e-02},
01085 {110.5,0.1588,148.3,-3.38,3.589e-02},
01086 {177.4,32.37,381.2,1.083,0.0}
01087 };
01088 static const double eth[7]={13.62,16.94,18.79,28.48,50.,110.5,538.};
01089 static const long l[7]={1,1,1,0,1,1,0};
01090
01091 DEBUG_ENTRY( "ofit()" );
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104 otot = 0.0;
01105
01106 for( int i=0; i < 7; i++ )
01107 {
01108 opart[i] = 0.0;
01109 }
01110
01111 for( int i=0; i < 7; i++ )
01112 {
01113 if( e >= eth[i] )
01114 {
01115 q = 5.5 - 0.5*y[i][3] + l[i];
01116
01117 x = e/y[i][0];
01118
01119 opart[i] = (float)(y[i][1]*(POW2(x - 1.0) + POW2(y[i][4]))/
01120 pow(x,q)/pow(1.0 + sqrt(x/y[i][2]),y[i][3]));
01121
01122 otot += opart[i];
01123 }
01124 }
01125
01126 DEBUG_EXIT( "ofit()" );
01127 return(otot);
01128 }
01129
01130
01131
01132
01133 static void OpacityCreate1Element(
01134
01135 long int nelem)
01136 {
01137 long int ihi,
01138 ip,
01139 ipop,
01140 limit,
01141 low,
01142 need,
01143 nelec,
01144 ion,
01145 nshell;
01146 double cs;
01147 double energy;
01148
01149 DEBUG_ENTRY( "OpacityCreate1Element()" );
01150
01151
01152 ASSERT( nelem >= 2 );
01153 ASSERT( nelem < LIMELM );
01154
01155
01156
01157 for( ion=0; ion < nelem; ion++ )
01158 {
01159
01160
01161 for( ip=0; ip < rfield.nupper; ip++ )
01162 {
01163 opac.opacity_abs[ip] = 0.;
01164 }
01165
01166
01167 nelec = nelem+1 - ion;
01168
01169
01170 for( nshell=0; nshell < Heavy.nsShells[nelem][ion]; nshell++ )
01171 {
01172
01173 opac.ipElement[nelem][ion][nshell][2] = opac.nOpacTot + 1;
01174
01175
01176 limit = opac.ipElement[nelem][ion][nshell][1];
01177
01178
01179 need = limit - opac.ipElement[nelem][ion][nshell][0] + 1;
01180
01181
01182 if( opac.nOpacTot + need > ndimOpacityStack )
01183 opacity_more_memory();
01184
01185
01186 low = opac.ipElement[nelem][ion][nshell][0];
01187 ihi = opac.ipElement[nelem][ion][nshell][1];
01188 ipop = opac.ipElement[nelem][ion][nshell][2];
01189
01190
01191
01192 ASSERT( low <= ihi || low<5 );
01193
01194
01195 for( ip=low-1; ip < ihi; ip++ )
01196 {
01197
01198 energy = MAX2(rfield.AnuOrg[ip]*EVRYD ,
01199 t_ADfA::Inst().ph1(nshell,nelec-1,nelem,0));
01200
01201
01202 cs = t_ADfA::Inst().phfit(nelem+1,nelec,nshell+1,energy);
01203
01204
01205
01206 opac.OpacStack[ip-low+ipop] = cs*1e-18;
01207
01208
01209 opac.opacity_abs[ip] += cs;
01210 }
01211
01212 opac.nOpacTot += ihi - low + 1;
01213
01214
01215 if( punch.lgPunPoint )
01216 {
01217 fprintf( punch.ipPoint, "%3ld%3ld%3ld%10.2e%10.2e%10.2e%10.2e\n",
01218 nelem, ion, nshell, rfield.anu[low-1], rfield.anu[ihi-1],
01219 opac.OpacStack[ipop-1], opac.OpacStack[ihi-low+ipop-1] );
01220 }
01221 }
01222
01223
01224 for(
01225 ip=opac.ipElement[nelem][ion][Heavy.nsShells[nelem][ion]-1][0]-1;
01226 ip < continuum.KshellLimit; ip++ )
01227 {
01228 ASSERT( opac.opacity_abs[ip] > 0. );
01229 }
01230
01231 }
01232
01233 DEBUG_EXIT( "OpacityCreate1Element()" );
01234 return;
01235 }
01236
01237
01238 static void opacity_more_memory(void)
01239 {
01240
01241 DEBUG_ENTRY( "opacity_more_memory()" );
01242
01243
01244 ndimOpacityStack *= 2;
01245 opac.OpacStack = (double *)REALLOC( opac.OpacStack , (size_t)ndimOpacityStack*sizeof(double) );
01246 fprintf( ioQQQ, " NOTE OpacityCreate1Element needed more opacity cells than ndimOpacityStack, please consider increasing it.\n" );
01247 fprintf( ioQQQ, " NOTE OpacityCreate1Element doubled memory allocation to %li.\n",ndimOpacityStack );
01248 lgRealloc = true;
01249
01250 DEBUG_EXIT( "opacity_more_memory()" );
01251 return;
01252 }
01253
01254
01255 static double Opacity_iso_photo_cs(
01256
01257 float energy ,
01258
01259 long ipISO ,
01260
01261 long nelem ,
01262
01263 long n )
01264 {
01265
01266 double thres;
01267 double crs=-DBL_MAX;
01268
01269 DEBUG_ENTRY( "Opacity_iso_photo_cs()" );
01270
01271 if( ipISO==ipH_LIKE )
01272 {
01273 double photon;
01274 if( n==0 )
01275 {
01276
01277
01278 thres = MAX2(energy*(float)EVRYD , t_ADfA::Inst().ph1(0,0,nelem,0));
01279 crs = t_ADfA::Inst().phfit(nelem+1,1,1,thres)* 1e-18;
01280
01281 ASSERT( crs > 0. && crs < 1e-10 );
01282 }
01283 else if( n==1 )
01284 {
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302 energy = MAX2( POW2( (float)(nelem+1) / 2.f ) , energy );
01303
01304
01305
01306 photon = energy / iso.xIsoLevNIonRyd[ipISO][nelem][n];
01307 photon = MAX2( photon , 1. + FLT_EPSILON*2. );
01308
01309
01310
01311 crs = H_photo_cs( photon , 2,0,nelem+1 );
01312
01313 ASSERT( crs > 0. && crs < 1e-10 );
01314 }
01315 else if( n==2 )
01316 {
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326 energy = MAX2( POW2( (float)(nelem+1) / 2.f ) , energy );
01327
01328
01329 photon = energy / iso.xIsoLevNIonRyd[ipISO][nelem][n];
01330 photon = MAX2( photon , 1. + FLT_EPSILON*2. );
01331
01332
01333
01334 crs = H_photo_cs( photon , 2,1,nelem+1 );
01335
01336 ASSERT( crs > 0. && crs < 1e-10 );
01337 }
01338 else
01339 {
01340
01341
01342
01343
01344
01345
01346 thres = MAX2( energy , iso.xIsoLevNIonRyd[ipISO][nelem][n]*1.001f );
01347
01348 crs = t_ADfA::Inst().hpfit(nelem+1,n,thres*EVRYD);
01349
01350 ASSERT( crs > 0. && crs < 1e-10 );
01351 }
01352 }
01353 else if( ipISO==ipHE_LIKE )
01354 {
01355 thres = MAX2( energy , iso.xIsoLevNIonRyd[ipISO][nelem][n]);
01356
01357 if( n >= iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] )
01358 {
01359 long int nup = iso.n_HighestResolved_max[ipHE_LIKE][nelem] + n + 1 -
01360 (iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem]);
01361
01362
01363
01364 crs = t_ADfA::Inst().hpfit(nelem,nup ,thres*EVRYD);
01365
01366 ASSERT(
01367 (energy < iso.xIsoLevNIonRyd[ipISO][nelem][n]*1.02) ||
01368 (crs > 0. && crs < 1e-10) );
01369 }
01370 else
01371 {
01372
01373
01374
01375
01376
01377
01378
01379 crs = He_cross_section( thres , n , nelem );
01380
01381 ASSERT( crs > 0. && crs < 1e-10 );
01382 }
01383 }
01384 else
01385 TotalInsanity();
01386
01387 DEBUG_EXIT( "Opacity_iso_photo_cs()" );
01388 return(crs);
01389 }
01390
01391
01392 static const int NCRS = 33;
01393
01394 static double hmiopc(double freq)
01395 {
01396 double energy,
01397 hmiopc_v,
01398 x,
01399 y;
01400 static double y2[NCRS];
01401 static double crs[NCRS]={0.,0.124,0.398,0.708,1.054,1.437,1.805,
01402 2.176,2.518,2.842,3.126,3.377,3.580,3.741,3.851,3.913,3.925,
01403 3.887,3.805,3.676,3.511,3.306,3.071,2.810,2.523,2.219,1.898,
01404 1.567,1.233,.912,.629,.39,.19};
01405 static double ener[NCRS]={0.,0.001459,0.003296,0.005256,0.007351,
01406 0.009595,0.01201,0.01460,0.01741,0.02044,0.02375,0.02735,0.03129,
01407 0.03563,0.04043,0.04576,0.05171,0.05841,0.06601,0.07469,0.08470,
01408 0.09638,0.1102,0.1268,0.1470,0.1723,0.2049,0.2483,0.3090,0.4001,
01409 0.5520,0.8557,1.7669};
01410 static bool lgFirst = true;
01411
01412 DEBUG_ENTRY( "hmiopc()" );
01413
01414
01415
01416
01417
01418
01419
01420 if( lgFirst )
01421 {
01422
01423 spline(ener,crs,NCRS,2e31,2e31,y2);
01424 lgFirst = false;
01425 }
01426
01427 energy = freq - 0.05552;
01428 if( energy < ener[0] || energy > ener[NCRS-1] )
01429 {
01430 hmiopc_v = 0.;
01431 }
01432 else
01433 {
01434 x = energy;
01435 splint(ener,crs,y2,NCRS,x,&y);
01436 hmiopc_v = y*1e-17;
01437 }
01438
01439 DEBUG_EXIT( "hmiopc()" );
01440 return( hmiopc_v );
01441 }
01442
01443
01444 static double rayleh(double ener)
01445 {
01446 double rayleh_v;
01447
01448 DEBUG_ENTRY( "rayleh()" );
01449
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461 if( ener < 0.05 )
01462 {
01463 rayleh_v = (8.41e-25*powi(ener,4) + 3.37e-24*powi(ener,6))*
01464 hydro.DampOnFac;
01465 }
01466
01467 else if( ener < 0.646 )
01468 {
01469 rayleh_v = (8.41e-25*powi(ener,4) + 3.37e-24*powi(ener,6) +
01470 4.71e-22*powi(ener,14))*hydro.DampOnFac;
01471 }
01472
01473 else if( ener >= 0.646 && ener < 1.0 )
01474 {
01475 rayleh_v = fabs(0.74959-ener);
01476 rayleh_v = 1.788e5/POW2(FR1RYD*MAX2(0.001,rayleh_v));
01477
01478 rayleh_v = MAX2(rayleh_v,1e-24)*hydro.DampOnFac;
01479 }
01480
01481 else
01482 {
01483 rayleh_v = 0.;
01484 }
01485
01486 DEBUG_EXIT( "rayleh()" );
01487 return( rayleh_v );
01488 }