00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "cddefines.h"
00011 #include "physconst.h"
00012 #include "lines_service.h"
00013 #include "iso.h"
00014 #include "secondaries.h"
00015 #include "taulines.h"
00016 #include "elementnames.h"
00017 #include "ionbal.h"
00018 #include "rt.h"
00019 #include "opacity.h"
00020 #include "yield.h"
00021 #include "dense.h"
00022 #include "he.h"
00023 #include "fe.h"
00024 #include "rfield.h"
00025 #include "oxy.h"
00026 #include "atomfeii.h"
00027 #include "atoms.h"
00028 #include "trace.h"
00029 #include "hmi.h"
00030 #include "heavy.h"
00031 #include "predcont.h"
00032 #include "atmdat.h"
00033 #include "ipoint.h"
00034 #include "h2.h"
00035 #include "continuum.h"
00036
00037
00038 STATIC long LimitSh(long int ion,
00039 long int nshell,
00040 long int nelem);
00041
00042 STATIC void ipShells(
00043
00044 long int nelem);
00045
00046
00047 STATIC void ContBandsCreate(
00048
00049
00050
00051 const char chFile[] );
00052
00053
00054 #define TwoS (1+ipISO)
00055
00056
00057 STATIC void fiddle(long int ipnt,
00058 double exact);
00059
00060 void ContCreatePointers(void)
00061 {
00062 char chLab[5];
00063 long int
00064 i,
00065 ion,
00066 ipHi,
00067 ipLo,
00068 ipISO,
00069 iWL_Ang,
00070 j,
00071 nelem,
00072 nshells;
00073 double energy,
00074 xnew;
00075
00076 static int nCalled = 0;
00077
00078 DEBUG_ENTRY( "ContCreatePointers()" );
00079
00080
00081
00082 iso_create();
00083
00084
00085 H2_Create();
00086
00087
00088
00089 if( nCalled > 0 )
00090 {
00091 if( trace.lgTrace )
00092 fprintf( ioQQQ, " ContCreatePointers called, not evaluating.\n" );
00093 return;
00094 }
00095 else
00096 {
00097 if( trace.lgTrace )
00098 fprintf( ioQQQ, " ContCreatePointers called first time.\n" );
00099 ++nCalled;
00100 }
00101
00102 for( i=0; i < rfield.nupper; i++ )
00103 {
00104
00105 strcpy( rfield.chContLabel[i], " ");
00106 strcpy( rfield.chLineLabel[i], " ");
00107 }
00108
00109
00110
00111
00112 for( nelem=0; nelem<LIMELM; ++nelem )
00113 {
00114 if( dense.lgElmtOn[nelem] )
00115 {
00116 for( ion=0; ion<LIMELM; ++ion )
00117 {
00118 for( nshells=0; nshells<7; ++nshells )
00119 {
00120 for( j=0; j<3; ++j )
00121 {
00122 opac.ipElement[nelem][ion][nshells][j] = INT_MIN;
00123 }
00124 }
00125 }
00126 }
00127 }
00128
00129
00130 opac.ipo3exc[0] = ipContEnergy(3.855,"O3ex");
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142 for( ipISO=ipH_LIKE; ipISO<=ipHE_LIKE; ++ipISO )
00143 {
00144
00145 for( nelem=ipISO; nelem < 2; nelem++ )
00146 {
00147
00148 sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO );
00149
00150
00151 iso.ipIsoLevNIonCon[ipISO][nelem][0] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][0],chLab);
00152 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00153 {
00154
00155 iso.ipIsoLevNIonCon[ipISO][nelem][ipHi] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][ipHi],chLab);
00156
00157
00158 for( ipLo=0; ipLo < ipHi; ipLo++ )
00159 {
00160 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
00161 continue;
00162
00163
00164
00165
00166
00167 if( Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD < continuum.filbnd[0] )
00168 continue;
00169
00170
00171 Transitions[ipISO][nelem][ipHi][ipLo].ipCont =
00172 ipLineEnergy(Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD , chLab ,
00173 iso.ipIsoLevNIonCon[ipISO][nelem][ipLo]);
00174 Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine =
00175 ipFineCont(Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD );
00176
00177 # ifndef NDEBUG
00178 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine > 0 )
00179 {
00180 realnum anuCoarse = rfield.anu[Transitions[ipISO][nelem][ipHi][ipLo].ipCont-1];
00181 realnum anuFine = rfield.fine_anu[Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine];
00182 realnum widCoarse = rfield.widflx[Transitions[ipISO][nelem][ipHi][ipLo].ipCont-1];
00183 realnum widFine = anuFine - rfield.fine_anu[Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine-1];
00184 realnum width = MAX2( widFine , widCoarse );
00185
00186
00187
00188
00189
00190 ASSERT( fabs(anuCoarse - anuFine) / anuCoarse <
00191 2.*width/anuCoarse);
00192 }
00193 # endif
00194 }
00195 }
00196 }
00197 }
00198
00199
00200
00201 nelem = ipHELIUM;
00202
00203 if( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] , "He 2" ) == 0)
00204 {
00205 strcpy( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]],
00206 rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] );
00207
00208 strcpy( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] , " ");
00209 }
00210 else if( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] , "H 1" ) == 0)
00211 {
00212
00213 strcpy( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]] , "He 2");
00214 }
00215 else
00216 {
00217 fprintf(ioQQQ," insanity heii pointer fix contcreatepointers\n");
00218 }
00219
00220 ++iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s];
00221 ++iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2p];
00222
00223
00224
00225
00226 secondaries.ipSecIon = ipoint(7.353);
00227
00228
00229
00230 continuum.KshellLimit = ipoint(continuum.EnergyKshell);
00231
00232
00233
00234 opac.ih2pnt[0] = ipContEnergy(0.502,"H2+ ");
00235 opac.ih2pnt[1] = ipoint(1.03);
00236
00237
00238
00239
00240
00241 i = ipContEnergy(0.0117, "PAH " );
00242
00243
00244 i = ipContEnergy(0.0147, "PAH " );
00245
00246
00247 i = ipContEnergy(0.028, "PAH " );
00248
00249
00250 i = ipContEnergy(0.0080, "PAH " );
00251
00252
00253 i = ipContEnergy(0.0077, "PAH " );
00254
00255
00256 i = ipContEnergy(0.0069, "PAH " );
00257
00258
00259
00260
00261
00262 he.ip374 = ipLineEnergy(1.92,"He 2",0);
00263 he.ip660 = ipLineEnergy(1.38,"He 2",0);
00264
00265
00266
00267 for( i=0; i < NPREDCONT; i++ )
00268 {
00269
00270 ipPredCont[i] = ipoint(EnrPredCont[i]) - 1;
00271 }
00272
00273
00274 rt.ipxry = ipoint(73.5);
00275
00276
00277 hmi.iphmin = ipContEnergy(0.05544,"H- ");
00278
00279
00280 opac.ipH2_photo_thresh = ipContEnergy( 15.4/EVRYD , "H2 ");
00281
00282 hmi.iheh1 = ipoint(1.6);
00283 hmi.iheh2 = ipoint(2.3);
00284
00285
00286 opac.ipCKshell = ipoint(20.6);
00287
00288
00289 opac.ippr = ipoint(7.51155e4) + 1;
00290
00291
00292 rfield.ipEnerGammaRay = ipoint(rfield.EnerGammaRay);
00293
00294 t_fe2ovr_la::Inst().init_pointers();
00295
00296
00297
00298
00299 fiddle(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1,0.99946);
00300 fiddle(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1,0.24986);
00301
00302 ASSERT( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1], "H 1" ) ==0 );
00303 ASSERT( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1], "H 1" ) ==0 );
00304
00305 fiddle(iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]-1,4.00000);
00306 ASSERT( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]-1], "He 2" ) ==0 );
00307
00308
00309 fiddle(opac.ipo3exc[0]-1,3.855);
00310
00311
00312 for( i=1; i<rfield.nupper-1; ++i )
00313 {
00314 rfield.widflx[i] = ((rfield.anu[i+1] - rfield.anu[i]) + (rfield.anu[i] - rfield.anu[i-1]))/2.f;
00315 }
00316
00317
00318
00319
00320 rfield.ipB_filter = ipoint( RYDLAM / WL_B_FILT );
00321
00322 rfield.ipV_filter = ipoint( RYDLAM / WL_V_FILT );
00323
00324
00325
00326 rfield.ipG0_TH85_lo = ipoint( 6.0 / EVRYD );
00327 rfield.ipG0_TH85_hi = ipoint( 13.6 / EVRYD );
00328
00329
00330 rfield.ipG0_DB96_lo = ipoint( RYDLAM / 1110. );
00331 rfield.ipG0_DB96_hi = ipoint( RYDLAM / 912. );
00332
00333
00334
00335 rfield.ipG0_spec_lo = ipoint( 6.0 / EVRYD );
00336 rfield.ipG0_spec_hi = ipoint( RYDLAM / 912. );
00337
00338
00339 rfield.ip1000A = ipoint( RYDLAM / 1000. );
00340
00341
00342 for( i=0; i < rfield.nupper; i++ )
00343 {
00344 rfield.AnuOrg[i] = rfield.anu[i];
00345 rfield.anusqr[i] = (realnum)sqrt(rfield.AnuOrg[i]);
00346 }
00347
00348
00349
00350
00351
00352
00353
00354 nelem = ipHYDROGEN;
00355 ion = 0;
00356
00357
00358 Heavy.nsShells[nelem][0] = 1;
00359
00360
00361 Heavy.ipHeavy[nelem][ion] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
00362 opac.ipElement[nelem][ion][0][0] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
00363
00364
00365 opac.ipElement[nelem][ion][0][1] = rfield.nupper;
00366
00367 if( dense.lgElmtOn[ipHELIUM] )
00368 {
00369
00370 nelem = ipHELIUM;
00371 ion = 0;
00372
00373
00374 Heavy.nsShells[nelem][0] = 1;
00375
00376
00377 Heavy.ipHeavy[nelem][ion] = iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0];
00378 opac.ipElement[nelem][ion][0][0] = iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0];
00379
00380
00381 opac.ipElement[nelem][ion][0][1] = rfield.nupper;
00382
00383
00384 ion = 1;
00385
00386 Heavy.nsShells[nelem][1] = 1;
00387
00388
00389 Heavy.ipHeavy[nelem][ion] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
00390 opac.ipElement[nelem][ion][0][0] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
00391
00392
00393 opac.ipElement[nelem][ion][0][1] = rfield.nupper;
00394 }
00395
00396
00397
00398 ASSERT( t_ADfA::Inst().ph1(2,5,5,0) > 0. );
00399
00400
00401
00402
00403 for( i=NISO; i<LIMELM; ++i )
00404 {
00405 if( dense.lgElmtOn[i])
00406 {
00407
00408 ipShells(i);
00409 }
00410 }
00411
00412
00413 Heavy.Valence_IP_Ryd[ipHYDROGEN][0] = t_ADfA::Inst().ph1(0,0,ipHYDROGEN,0)/EVRYD* 0.9998787;
00414 Heavy.Valence_IP_Ryd[ipHELIUM][0] = t_ADfA::Inst().ph1(0,1,ipHELIUM,0)/EVRYD* 0.9998787;
00415 Heavy.Valence_IP_Ryd[ipHELIUM][1] = t_ADfA::Inst().ph1(0,0,ipHELIUM,0)/EVRYD* 0.9998787;
00416 for( nelem=2; nelem<LIMELM; ++nelem )
00417 {
00418 Heavy.Valence_IP_Ryd[nelem][nelem-1] = t_ADfA::Inst().ph1(0,1,nelem,0)/EVRYD* 0.9998787;
00419 Heavy.Valence_IP_Ryd[nelem][nelem] = t_ADfA::Inst().ph1(0,0,nelem,0)/EVRYD* 0.9998787;
00420 if( dense.lgElmtOn[nelem])
00421 {
00422
00423 for( j=0; j<=nelem; ++j )
00424 {
00425 ASSERT( Heavy.Valence_IP_Ryd[nelem][j] > 0.05 );
00426 }
00427 for( j=0; j<nelem; ++j )
00428 {
00429 ASSERT( Heavy.Valence_IP_Ryd[nelem][j] < Heavy.Valence_IP_Ryd[nelem][j+1]);
00430 }
00431 }
00432 }
00433
00434
00435 for( nelem=0; nelem<LIMELM; ++nelem )
00436 {
00437 if( dense.lgElmtOn[nelem])
00438 {
00439 for( ion=0; ion<nelem+1; ++ion )
00440 {
00441
00442 energy = sqrt( Heavy.Valence_IP_Ryd[nelem][ion] * EN1RYD * ELECTRON_MASS * SPEEDLIGHT * SPEEDLIGHT ) / EN1RYD;
00443
00444 ionbal.ipCompRecoil[nelem][ion] = ipoint( energy );
00445 }
00446 }
00447 }
00448
00449
00450
00451 oxy.i2d = ipoint(1.242);
00452 oxy.i2p = ipoint(1.367);
00453 opac.ipo1exc[0] = ipContEnergy(0.856,"O1ex");
00454 opac.ipo1exc[1] = ipoint(2.0);
00455
00456
00457
00458 opac.ipo3exc[1] = ipoint(5.0);
00459
00460
00461 opac.ipo3exc3[0] = ipContEnergy(3.646,"O3ex");
00462 opac.ipo3exc3[1] = ipoint(5.0);
00463
00464
00465
00470 atoms.ipoiex[0] = ipoint(0.7005);
00471 atoms.ipoiex[1] = ipoint(0.10791);
00472 atoms.ipoiex[2] = ipoint(0.06925);
00473 atoms.ipoiex[3] = ipoint(0.01151);
00474 atoms.ipoiex[4] = ipoint(0.01999);
00475
00476
00477
00478
00479
00480 opac.in1[0] = ipContEnergy(0.893,"N1ex");
00481
00482
00483 opac.in1[1] = ipoint(2.);
00484
00485 if( (trace.lgTrace && trace.lgConBug) || (trace.lgTrace && trace.lgPointBug) )
00486 {
00487 fprintf( ioQQQ, " ContCreatePointers:%ld energy cells used. N(1R):%4ld N(1.8):%4ld N(4Ryd):%4ld N(O3)%4ld N(x-ray):%5ld N(rcoil)%5ld\n",
00488 rfield.nupper,
00489 iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s],
00490 iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][ipH1s],
00491 iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s],
00492 opac.ipo3exc[0],
00493 opac.ipCKshell,
00494 ionbal.ipCompRecoil[ipHYDROGEN][0] );
00495
00496
00497 fprintf( ioQQQ, " ContCreatePointers: ipEnerGammaRay: %5ld IPPRpari produc%5ld\n",
00498 rfield.ipEnerGammaRay, opac.ippr );
00499
00500 fprintf( ioQQQ, " ContCreatePointers: H pointers;" );
00501 for( i=0; i <= 6; i++ )
00502 {
00503 fprintf( ioQQQ, "%4ld%4ld", i, iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][i] );
00504 }
00505 fprintf( ioQQQ, "\n" );
00506
00507 fprintf( ioQQQ, " ContCreatePointers: Oxy pnters;" );
00508
00509 for( i=1; i <= 8; i++ )
00510 {
00511 fprintf( ioQQQ, "%4ld%4ld", i, Heavy.ipHeavy[ipOXYGEN][i-1] );
00512 }
00513 fprintf( ioQQQ, "\n" );
00514 }
00515
00516
00517
00518 opac.ipmgex = ipoint(0.779);
00519
00520
00521 opac.ica2ex[0] = ipContEnergy(0.72,"Ca2x");
00522 opac.ica2ex[1] = ipoint(1.);
00523
00524
00525 fe.ipfe10 = ipoint(2.605);
00526
00527
00528 fe.pfe10 = (realnum)(2.00e-18/rfield.widflx[fe.ipfe10-1]);
00529
00530
00531 fe.pfe11a = (realnum)(4.54e-19/rfield.widflx[fe.ipfe10-1]);
00532
00533
00534 fe.pfe11b = (realnum)(2.75e-19/rfield.widflx[fe.ipfe10-1]);
00535 fe.pfe14 = (realnum)(1.15e-18/rfield.widflx[fe.ipfe10-1]);
00536
00537
00538
00539
00540
00541 for( i=1; i <= nLevel1; i++ )
00542 {
00543
00544 chIonLbl(chLab, &TauLines[i]);
00545
00546 TauLines[i].ipCont =
00547 ipLineEnergy(TauLines[i].EnergyWN * WAVNRYD, chLab ,0);
00548 TauLines[i].Emis->ipFine =
00549 ipFineCont(TauLines[i].EnergyWN * WAVNRYD );
00550
00551
00552
00553
00554
00555 if( TauLines[i].Emis->gf > 0. )
00556 {
00557
00558
00559 TauLines[i].Emis->Aul =
00560 (realnum)(eina(TauLines[i].Emis->gf,
00561 TauLines[i].EnergyWN,TauLines[i].Hi->g));
00562 }
00563 else if( TauLines[i].Emis->Aul > 0. )
00564 {
00565
00566
00567 TauLines[i].Emis->gf =
00568 (realnum)(GetGF(TauLines[i].Emis->Aul,
00569 TauLines[i].EnergyWN,TauLines[i].Hi->g));
00570 }
00571 else
00572 {
00573 fprintf( ioQQQ, " level 1 line does not have valid gf or A\n" );
00574 fprintf( ioQQQ, " This is ContCreatePointers\n" );
00575 cdEXIT(EXIT_FAILURE);
00576 }
00577
00578
00579 TauLines[i].Emis->dampXvel =
00580 (realnum)(TauLines[i].Emis->Aul / TauLines[i].EnergyWN/PI4);
00581
00582
00583 TauLines[i].Emis->opacity =
00584 (realnum)(abscf(TauLines[i].Emis->gf,
00585 TauLines[i].EnergyWN,TauLines[i].Lo->g));
00586
00587
00588
00589
00590 TauLines[i].EnergyK =
00591 (realnum)(T1CM)*TauLines[i].EnergyWN;
00592
00593
00594 TauLines[i].EnergyErg =
00595 (realnum)(ERG1CM)*TauLines[i].EnergyWN;
00596
00597
00598 ASSERT( TauLines[i].WLAng > 0 );
00599 }
00600
00601
00602 for( i=0; i <linesAdded2; i++)
00603 {
00604 atmolEmis[i].gf = (realnum)GetGF(atmolEmis[i].Aul,
00605 atmolEmis[i].tran->EnergyWN,atmolEmis[i].tran->Hi->g);
00606 atmolEmis[i].dampXvel = (realnum)(atmolEmis[i].Aul/
00607 atmolEmis[i].tran->EnergyWN/PI4);
00608 atmolEmis[i].damp = -1000.0;
00609
00610 strncpy(chLab,atmolEmis[i].tran->Hi->chLabel,4);
00611 chLab[4]='\0';
00612
00613 atmolEmis[i].tran->ipCont =
00614 ipLineEnergy(atmolEmis[i].tran->EnergyWN * WAVNRYD, chLab ,0);
00615 atmolEmis[i].ipFine = ipFineCont(atmolEmis[i].tran->EnergyWN * WAVNRYD );
00616
00617 atmolEmis[i].opacity =
00618 (realnum)(abscf(atmolEmis[i].gf,atmolEmis[i].tran->EnergyWN,
00619 atmolEmis[i].tran->Lo->g));
00620
00621 atmolEmis[i].tran->EnergyK =
00622 (realnum)(T1CM)*atmolEmis[i].tran->EnergyWN;
00623
00624 atmolEmis[i].tran->EnergyErg =
00625 (realnum)(ERG1CM)*atmolEmis[i].tran->EnergyWN ;
00626 }
00627
00628
00629
00630 H2_ContPoint();
00631
00632 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00633 {
00634
00635 for( nelem=2; nelem < LIMELM; nelem++ )
00636 {
00637 if( dense.lgElmtOn[nelem])
00638 {
00639
00640 sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO );
00641
00642
00643
00644 iso.ipIsoLevNIonCon[ipISO][nelem][0] =
00645 ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][0],chLab);
00646
00647 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00648 {
00649
00650 iso.ipIsoLevNIonCon[ipISO][nelem][ipHi] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][ipHi],chLab);
00651
00652
00653 for( ipLo=0; ipLo < ipHi; ipLo++ )
00654 {
00655
00656 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
00657 continue;
00658
00659
00660
00661
00662
00663 if( Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD < continuum.filbnd[0] )
00664 continue;
00665
00666 Transitions[ipISO][nelem][ipHi][ipLo].ipCont =
00667 ipLineEnergy(Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD , chLab ,
00668 iso.ipIsoLevNIonCon[ipISO][nelem][ipLo]);
00669 Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine =
00670 ipFineCont(Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD );
00671 }
00672 }
00673 iso.ipIsoLevNIonCon[ipISO][nelem][0] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][0],chLab);
00674 }
00675 }
00676 }
00677 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00678 {
00679
00680 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00681 {
00682 if( dense.lgElmtOn[nelem])
00683 {
00684 ipLo = 0;
00685
00686 for( ipHi=2; ipHi < iso.nLyman[ipISO]; ipHi++ )
00687 {
00688
00689
00690 ExtraLymanLines[ipISO][nelem][ipHi].ipCont =
00691 ipLineEnergy(ExtraLymanLines[ipISO][nelem][ipHi].EnergyWN * WAVNRYD , chLab ,
00692 iso.ipIsoLevNIonCon[ipISO][nelem][ipLo]);
00693
00694 ExtraLymanLines[ipISO][nelem][ipHi].Emis->ipFine =
00695 ipFineCont(ExtraLymanLines[ipISO][nelem][ipHi].EnergyWN * WAVNRYD );
00696
00697 }
00698
00699 if( iso.lgDielRecom[ipISO] )
00700 {
00701 ASSERT( ipISO>ipH_LIKE );
00702 for( ipHi=0; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
00703 {
00704
00705 SatelliteLines[ipISO][nelem][ipHi].ipCont = ipLineEnergy(
00706 SatelliteLines[ipISO][nelem][ipHi].EnergyWN * WAVNRYD , chLab ,
00707 0);
00708
00709 SatelliteLines[ipISO][nelem][ipHi].Emis->ipFine =
00710 ipFineCont(SatelliteLines[ipISO][nelem][ipHi].EnergyWN * WAVNRYD );
00711 }
00712 }
00713 }
00714 }
00715 }
00716
00717
00718
00719 ipISO = ipHYDROGEN;
00720
00721 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00722 {
00723 if( dense.lgElmtOn[nelem])
00724 {
00725
00726 Transitions[ipISO][nelem][ipH2p][ipH2s].ipCont = -1;
00727
00728 Transitions[ipISO][nelem][ipH2s][ipH1s].ipCont = -1;
00729
00730 }
00731 }
00732
00733 fixit();
00734
00735 ipISO = ipHELIUM;
00736
00737 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00738 {
00739 if( dense.lgElmtOn[nelem])
00740 {
00741
00742 Transitions[ipISO][nelem][ipHe2s1S][ipHe1s1S].ipCont = -1;
00743 Transitions[ipISO][nelem][ipHe2s1S][ipHe1s1S].Emis->ipFine = -1;
00744
00745 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00746 {
00747 for( ipLo=0; ipLo < ipHi; ipLo++ )
00748 {
00749 if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 )
00750 continue;
00751
00752 if( fabs(Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul - iso.SmallA) < SMALLFLOAT )
00753 {
00754
00755 Transitions[ipISO][nelem][ipHi][ipLo].ipCont = -1;
00756 Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine = -1;
00757 }
00758 }
00759 }
00760 }
00761 }
00762
00763
00764 for( i=0; i < nCORotate; i++ )
00765 {
00766
00767 C12O16Rotate[i].EnergyK =
00768 (realnum)(T1CM)*C12O16Rotate[i].EnergyWN;
00769 C13O16Rotate[i].EnergyK =
00770 (realnum)(T1CM)*C13O16Rotate[i].EnergyWN;
00771
00772
00773 C12O16Rotate[i].EnergyErg =
00774 (realnum)(ERG1CM)*C12O16Rotate[i].EnergyWN;
00775 C13O16Rotate[i].EnergyErg =
00776 (realnum)(ERG1CM)*C13O16Rotate[i].EnergyWN;
00777
00778
00779 chIonLbl(chLab, &C12O16Rotate[i]);
00780 chIonLbl(chLab, &C13O16Rotate[i]);
00781
00782 C12O16Rotate[i].ipCont =
00783 ipLineEnergy(C12O16Rotate[i].EnergyWN * WAVNRYD, "12CO" ,0);
00784 C12O16Rotate[i].Emis->ipFine =
00785 ipFineCont(C12O16Rotate[i].EnergyWN * WAVNRYD );
00786 C13O16Rotate[i].ipCont =
00787 ipLineEnergy(C13O16Rotate[i].EnergyWN * WAVNRYD, "13CO" ,0);
00788 C13O16Rotate[i].Emis->ipFine =
00789 ipFineCont(C13O16Rotate[i].EnergyWN * WAVNRYD );
00790
00791 if( C12O16Rotate[i].Emis->gf > 0. )
00792 {
00793
00794
00795 C12O16Rotate[i].Emis->Aul =
00796 (realnum)(eina(C12O16Rotate[i].Emis->gf,
00797 C12O16Rotate[i].EnergyWN,C12O16Rotate[i].Hi->g));
00798 }
00799 else if( C12O16Rotate[i].Emis->Aul > 0. )
00800 {
00801
00802
00803 C12O16Rotate[i].Emis->gf =
00804 (realnum)(GetGF(C12O16Rotate[i].Emis->Aul,
00805 C12O16Rotate[i].EnergyWN,C12O16Rotate[i].Hi->g));
00806 }
00807 else
00808 {
00809 fprintf( ioQQQ, " 12CO line does not have valid gf or A\n" );
00810 fprintf( ioQQQ, " This is ContCreatePointers\n" );
00811 cdEXIT(EXIT_FAILURE);
00812 }
00813 if( C13O16Rotate[i].Emis->gf > 0. )
00814 {
00815
00816
00817 C13O16Rotate[i].Emis->Aul =
00818 (realnum)(eina(C13O16Rotate[i].Emis->gf,
00819 C13O16Rotate[i].EnergyWN,C13O16Rotate[i].Hi->g));
00820 }
00821 else if( C13O16Rotate[i].Emis->Aul > 0. )
00822 {
00823
00824
00825 C13O16Rotate[i].Emis->gf =
00826 (realnum)(GetGF(C13O16Rotate[i].Emis->Aul,
00827 C13O16Rotate[i].EnergyWN,C13O16Rotate[i].Hi->g));
00828 }
00829 else
00830 {
00831 fprintf( ioQQQ, " 13CO line does not have valid gf or A\n" );
00832 fprintf( ioQQQ, " This is ContCreatePointers\n" );
00833 cdEXIT(EXIT_FAILURE);
00834 }
00835
00836
00837 ASSERT( C12O16Rotate[i].WLAng > 0 );
00838 ASSERT( C13O16Rotate[i].WLAng > 0 );
00839
00840 C12O16Rotate[i].Emis->dampXvel =
00841 (realnum)(C12O16Rotate[i].Emis->Aul/
00842 C12O16Rotate[i].EnergyWN/PI4);
00843
00844 C13O16Rotate[i].Emis->dampXvel =
00845 (realnum)(C13O16Rotate[i].Emis->Aul/
00846 C13O16Rotate[i].EnergyWN/PI4);
00847
00848
00849 C12O16Rotate[i].Emis->opacity =
00850 (realnum)(abscf(C12O16Rotate[i].Emis->gf,
00851 C12O16Rotate[i].EnergyWN,C12O16Rotate[i].Lo->g));
00852 C13O16Rotate[i].Emis->opacity =
00853 (realnum)(abscf(C13O16Rotate[i].Emis->gf,
00854 C13O16Rotate[i].EnergyWN,C13O16Rotate[i].Lo->g));
00855 }
00856
00857
00858 for( i=0; i<nUTA; ++i )
00859 {
00860 if( UTALines[i].Emis->Aul > 0. )
00861 {
00862 UTALines[i].Emis->dampXvel =
00863 (realnum)(UTALines[i].Emis->Aul/ UTALines[i].EnergyWN/PI4);
00864
00865
00866 UTALines[i].Emis->opacity =
00867 (realnum)(abscf( UTALines[i].Emis->gf, UTALines[i].EnergyWN, UTALines[i].Lo->g));
00868
00869
00870 UTALines[i].EnergyK =
00871 (realnum)(T1CM*UTALines[i].EnergyWN);
00872
00873
00874 UTALines[i].EnergyErg =
00875 (realnum)(ERG1CM*UTALines[i].EnergyWN);
00876
00877
00878 chIonLbl(chLab, &UTALines[i]);
00879
00880
00881 UTALines[i].ipCont = ipLineEnergy(UTALines[i].EnergyWN * WAVNRYD , chLab,0 );
00882 UTALines[i].Emis->ipFine = ipFineCont(UTALines[i].EnergyWN * WAVNRYD );
00883
00884
00885
00886
00887 double thresh = Heavy.Valence_IP_Ryd[UTALines[i].Hi->nelem-1][UTALines[i].Hi->IonStg-1] *EN1RYD;
00888 UTALines[i].Coll.heat = (UTALines[i].EnergyErg-thresh);
00889 ASSERT( UTALines[i].Coll.heat> 0. );
00890 }
00891 }
00892
00893
00894 for( i=0; i < nWindLine; i++ )
00895 {
00896
00897
00898 TauLine2[i].Emis->Aul =
00899 (realnum)(eina(TauLine2[i].Emis->gf,
00900 TauLine2[i].EnergyWN,TauLine2[i].Hi->g));
00901
00902
00903 TauLine2[i].Emis->dampXvel =
00904 (realnum)(TauLine2[i].Emis->Aul/
00905 TauLine2[i].EnergyWN/PI4);
00906
00907
00908 TauLine2[i].Emis->opacity =
00909 (realnum)(abscf(TauLine2[i].Emis->gf,
00910 TauLine2[i].EnergyWN,TauLine2[i].Lo->g));
00911
00912
00913 TauLine2[i].EnergyK =
00914 (realnum)(T1CM*TauLine2[i].EnergyWN);
00915
00916
00917 TauLine2[i].EnergyErg =
00918 (realnum)(ERG1CM*TauLine2[i].EnergyWN);
00919
00920
00921 chIonLbl(chLab, &TauLine2[i]);
00922
00923
00924 TauLine2[i].ipCont = ipLineEnergy(TauLine2[i].EnergyWN * WAVNRYD , chLab,0 );
00925 TauLine2[i].Emis->ipFine = ipFineCont(TauLine2[i].EnergyWN * WAVNRYD );
00926
00927
00928 }
00929
00930
00931 for( i=0; i < nHFLines; i++ )
00932 {
00933
00934
00935
00936
00937
00938
00939 HFLines[i].Emis->dampXvel =
00940 (realnum)(HFLines[i].Emis->Aul/
00941 HFLines[i].EnergyWN/PI4);
00942 HFLines[i].Emis->damp = 1e-20f;
00943
00944
00945 HFLines[i].Emis->opacity =
00946 (realnum)(abscf(HFLines[i].Emis->gf,
00947 HFLines[i].EnergyWN,
00948 HFLines[i].Lo->g));
00949
00950
00951
00952
00953
00954 HFLines[i].EnergyK =
00955 (realnum)(T1CM*HFLines[i].EnergyWN);
00956
00957
00958 HFLines[i].EnergyErg =
00959 (realnum)(ERG1CM*HFLines[i].EnergyWN);
00960
00961
00962 chIonLbl(chLab, &HFLines[i]);
00963
00964
00965 HFLines[i].ipCont = ipLineEnergy(HFLines[i].EnergyWN * WAVNRYD , chLab,0 );
00966 HFLines[i].Emis->ipFine = ipFineCont(HFLines[i].EnergyWN * WAVNRYD );
00967 }
00968
00969
00970
00971 FeIIPoint();
00972
00973
00974 for( i=0; i < t_yield::Inst().nlines(); ++i )
00975 {
00976 strcpy( chLab , elementnames.chElementSym[t_yield::Inst().nelem(i)] );
00977 strcat( chLab , elementnames.chIonStage[t_yield::Inst().ion_emit(i)] );
00978
00979 t_yield::Inst().set_ipoint( i, ipLineEnergy( t_yield::Inst().energy(i) , chLab , 0 ) );
00980 }
00981
00982
00983
00984
00985
00986
00987
00988
00989 iso.ipSym2nu.reserve( NISO );
00990
00991
00992 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00993 {
00994 iso.ipSym2nu.reserve( ipISO, LIMELM );
00995
00996
00997 for( nelem=ipISO; nelem<LIMELM; ++nelem )
00998 {
00999 if( dense.lgElmtOn[nelem] )
01000 {
01001 double E2nu = Transitions[ipISO][nelem][TwoS][0].EnergyWN * WAVNRYD;
01002
01003
01004 iso.ipTwoPhoE[ipISO][nelem] = ipoint(E2nu);
01005
01006 iso.ipHalfTwoPhoE[ipISO][nelem] = ipoint(E2nu / 2.);
01007 while( rfield.anu[iso.ipTwoPhoE[ipISO][nelem]] > E2nu )
01008 {
01009 --iso.ipTwoPhoE[ipISO][nelem];
01010 }
01011 while( rfield.anu[iso.ipHalfTwoPhoE[ipISO][nelem]] > E2nu / 2. )
01012 {
01013 --iso.ipHalfTwoPhoE[ipISO][nelem];
01014 }
01015
01016 iso.ipSym2nu.reserve( ipISO, nelem, iso.ipTwoPhoE[ipISO][nelem] );
01017 }
01018 }
01019 }
01020
01021 iso.ipSym2nu.alloc();
01022 iso.As2nu.alloc( iso.ipSym2nu.clone() );
01023
01024
01025 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
01026 {
01027
01028 for( nelem=ipISO; nelem<LIMELM; ++nelem )
01029 {
01030 if( dense.lgElmtOn[nelem] )
01031 {
01032 double E2nu = Transitions[ipISO][nelem][TwoS][0].EnergyWN * WAVNRYD;
01033 double Aul = Transitions[ipISO][nelem][TwoS][0].Emis->Aul;
01034 double SumShapeFunction = 0., Renorm= 0.;
01035
01036
01037 for( i=0; i < iso.ipTwoPhoE[ipISO][nelem]; i++ )
01038 {
01039
01040 energy = E2nu - rfield.anu[i];
01041
01042
01043 energy = MAX2( energy, rfield.anu[0] + rfield.widflx[0]/2. );
01044
01045 iso.ipSym2nu[ipISO][nelem][i] = ipoint(energy);
01046 while( rfield.anu[iso.ipSym2nu[ipISO][nelem][i]] > E2nu ||
01047 iso.ipSym2nu[ipISO][nelem][i] >= iso.ipTwoPhoE[ipISO][nelem])
01048 {
01049 --iso.ipSym2nu[ipISO][nelem][i];
01050 }
01051 ASSERT( iso.ipSym2nu[ipISO][nelem][i] >= 0 );
01052 }
01053
01054
01055
01056 for( i=0; i<iso.ipTwoPhoE[ipISO][nelem]; i++ )
01057 {
01058 double ShapeFunction;
01059
01060 ASSERT( rfield.anu[i]<=E2nu );
01061
01062 ShapeFunction = atmdat_2phot_shapefunction( rfield.AnuOrg[i]/E2nu, ipISO, nelem )*rfield.widflx[i]/E2nu;
01063 SumShapeFunction += ShapeFunction;
01064
01065
01066
01067 iso.As2nu[ipISO][nelem][i] = (realnum)( Aul * ShapeFunction );
01068 }
01069
01070
01071
01072 Renorm = 1./SumShapeFunction;
01073
01074 for( i=0; i<iso.ipTwoPhoE[ipISO][nelem]; i++ )
01075 {
01076 iso.As2nu[ipISO][nelem][i] *= (realnum)Renorm;
01077 }
01078
01079
01080 ASSERT( fabs( SumShapeFunction*Renorm - 1. ) < 0.00001 );
01081 }
01082 }
01083 }
01084
01085 {
01086
01087 enum {DEBUG_LOC=false};
01088 if( DEBUG_LOC )
01089 {
01090 # define NCRS 21
01091 double limit,ener[NCRS]={
01092 0., 0.03738, 0.07506, 0.1124, 0.1498, 0.1875,
01093 0.225, 0.263, 0.300, 0.3373, 0.375, 0.4127,
01094 0.4500, 0.487, 0.525, 0.5625, 0.6002, 0.6376,
01095 0.6749, 0.7126, 0.75};
01096
01097 nelem = ipHYDROGEN;
01098 ipISO = ipHYDROGEN;
01099
01100 limit = iso.ipTwoPhoE[ipISO][nelem];
01101
01102 for( i=0; i < NCRS; i++ )
01103 {
01104 fprintf(ioQQQ,"%.3e\t%.3e\n", ener[i] ,
01105 atmdat_2phot_shapefunction( ener[i]/0.75, ipISO, nelem ) );
01106 }
01107
01108 xnew = 0.;
01110 for( i=0; i < limit; i++ )
01111 {
01112 double fach = iso.As2nu[ipISO][nelem][i]/2.*rfield.anu2[i]/rfield.widflx[i]*EN1RYD;
01113 fprintf(ioQQQ,"%.3e\t%.3e\t%.3e\n",
01114 rfield.anu[i] ,
01115 iso.As2nu[ipISO][nelem][i] / rfield.widflx[i] ,
01116 fach );
01117 xnew += iso.As2nu[ipISO][nelem][i];
01118 }
01119 fprintf(ioQQQ," sum is %.3e\n", xnew );
01120 cdEXIT(EXIT_FAILURE);
01121 }
01122 }
01123
01124 {
01125
01126 enum {DEBUG_LOC=false};
01127 if( DEBUG_LOC )
01128 {
01129 for( i=0; i<11; ++i )
01130 {
01131 char chLsav[10];
01132 TauDummy.WLAng = (realnum)(PI * pow(10.,(double)i));
01133 strcpy( chLsav, chLineLbl(&TauDummy) );
01134 fprintf(ioQQQ,"%.2f\t%s\n", TauDummy.WLAng , chLsav );
01135 }
01136 cdEXIT(EXIT_FAILURE);
01137 }
01138 }
01139
01140
01141 if( trace.lgTrLine )
01142 {
01143 fprintf( ioQQQ, " WL(Ang) E(RYD) IP gl gu gf A damp abs K\n" );
01144 for( i=1; i <= nLevel1; i++ )
01145 {
01146 strcpy( chLab, chLineLbl(&TauLines[i]) );
01147 iWL_Ang = (long)TauLines[i].WLAng;
01148 if( iWL_Ang > 1000000 )
01149 {
01150 iWL_Ang /= 10000;
01151 }
01152 else if( iWL_Ang > 10000 )
01153 {
01154 iWL_Ang /= 1000;
01155 }
01156
01157 fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
01158 chLab, iWL_Ang, RYDLAM/TauLines[i].WLAng,
01159 TauLines[i].ipCont, (long)(TauLines[i].Lo->g),
01160 (long)(TauLines[i].Hi->g), TauLines[i].Emis->gf,
01161 TauLines[i].Emis->Aul, TauLines[i].Emis->dampXvel,
01162 TauLines[i].Emis->opacity );
01163 }
01164
01165
01166 for( i=0; i <linesAdded2; i++)
01167 {
01168 strcpy( chLab, chLineLbl(atmolEmis[i].tran) );
01169
01170 iWL_Ang = (long)atmolEmis[i].tran->WLAng;
01171
01172 if( iWL_Ang > 1000000 )
01173 {
01174 iWL_Ang /= 10000;
01175 }
01176 else if( iWL_Ang > 10000 )
01177 {
01178 iWL_Ang /= 1000;
01179 }
01180 fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
01181 chLab, iWL_Ang, RYDLAM/atmolEmis[i].tran->WLAng,
01182 atmolEmis[i].tran->ipCont, (long)(atmolEmis[i].tran->Lo->g),
01183 (long)(atmolEmis[i].tran->Hi->g),atmolEmis[i].gf,
01184 atmolEmis[i].Aul,atmolEmis[i].dampXvel,
01185 atmolEmis[i].opacity);
01186 }
01187
01188
01189 for( i=0; i < nCORotate; i++ )
01190 {
01191 strcpy( chLab, chLineLbl(&C12O16Rotate[i]) );
01192
01193 iWL_Ang = (long)C12O16Rotate[i].WLAng;
01194
01195 if( iWL_Ang > 1000000 )
01196 {
01197 iWL_Ang /= 10000;
01198 }
01199 else if( iWL_Ang > 10000 )
01200 {
01201 iWL_Ang /= 1000;
01202 }
01203 fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
01204 chLab, iWL_Ang, RYDLAM/C12O16Rotate[i].WLAng,
01205 C12O16Rotate[i].ipCont, (long)(C12O16Rotate[i].Lo->g),
01206 (long)(C12O16Rotate[i].Hi->g), C12O16Rotate[i].Emis->gf,
01207 C12O16Rotate[i].Emis->Aul, C12O16Rotate[i].Emis->dampXvel,
01208 C12O16Rotate[i].Emis->opacity );
01209 }
01210
01211
01212 for( i=0; i < nCORotate; i++ )
01213 {
01214 strcpy( chLab, chLineLbl(&C13O16Rotate[i]) );
01215
01216 iWL_Ang = (long)C13O16Rotate[i].WLAng;
01217
01218 if( iWL_Ang > 1000000 )
01219 {
01220 iWL_Ang /= 10000;
01221 }
01222 else if( iWL_Ang > 10000 )
01223 {
01224 iWL_Ang /= 1000;
01225 }
01226 fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
01227 chLab, iWL_Ang, RYDLAM/C13O16Rotate[i].WLAng,
01228 C13O16Rotate[i].ipCont, (long)(C13O16Rotate[i].Lo->g),
01229 (long)(C13O16Rotate[i].Hi->g), C13O16Rotate[i].Emis->gf,
01230 C13O16Rotate[i].Emis->Aul, C13O16Rotate[i].Emis->dampXvel,
01231 C13O16Rotate[i].Emis->opacity );
01232 }
01233
01234 for( i=0; i < nWindLine; i++ )
01235 {
01236 strcpy( chLab, chLineLbl(&TauLine2[i]) );
01237
01238 iWL_Ang = (long)TauLine2[i].WLAng;
01239
01240 if( iWL_Ang > 1000000 )
01241 {
01242 iWL_Ang /= 10000;
01243 }
01244 else if( iWL_Ang > 10000 )
01245 {
01246 iWL_Ang /= 1000;
01247 }
01248 fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
01249 chLab, iWL_Ang, RYDLAM/TauLine2[i].WLAng,
01250 TauLine2[i].ipCont, (long)(TauLine2[i].Lo->g),
01251 (long)(TauLine2[i].Hi->g), TauLine2[i].Emis->gf,
01252 TauLine2[i].Emis->Aul, TauLine2[i].Emis->dampXvel,
01253 TauLine2[i].Emis->opacity );
01254 }
01255 for( i=0; i < nHFLines; i++ )
01256 {
01257 strcpy( chLab, chLineLbl(&HFLines[i]) );
01258
01259 iWL_Ang = (long)HFLines[i].WLAng;
01260
01261 if( iWL_Ang > 1000000 )
01262 {
01263 iWL_Ang /= 10000;
01264 }
01265 else if( iWL_Ang > 10000 )
01266 {
01267 iWL_Ang /= 1000;
01268 }
01269 fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
01270 chLab, iWL_Ang, RYDLAM/HFLines[i].WLAng,
01271 HFLines[i].ipCont, (long)(HFLines[i].Lo->g),
01272 (long)(HFLines[i].Hi->g), HFLines[i].Emis->gf,
01273 HFLines[i].Emis->Aul, HFLines[i].Emis->dampXvel,
01274 HFLines[i].Emis->opacity );
01275 }
01276 }
01277
01278
01279 if( !rt.lgFstOn )
01280 {
01281 for( i=1; i <= nLevel1; i++ )
01282 {
01283 if( TauLines[i].EnergyWN < 10000. )
01284 {
01285 TauLines[i].Emis->opacity = 0.;
01286 }
01287 }
01288
01289
01290 for( i=0; i <linesAdded2; i++)
01291 {
01292 if(atmolEmis[i].tran->EnergyWN < 10000. )
01293 {
01294 atmolEmis[i].opacity = 0.;
01295 }
01296 }
01297 }
01298
01299
01300 ContBandsCreate( "" );
01301
01302
01303
01304 lgLinesAdded = true;
01305 lgStatesAdded = true;
01306
01307 return;
01308 }
01309
01310
01311
01312
01313 STATIC void fiddle(long int ipnt,
01314 double exact)
01315 {
01316 realnum Ehi,
01317 Elo,
01318 OldEner;
01319
01320 DEBUG_ENTRY( "fiddle()" );
01321
01322
01323 ASSERT( ipnt >= 0 );
01324 ASSERT( ipnt < rfield.nupper-1 );
01325
01326
01327 Ehi = rfield.anu[ipnt] + rfield.widflx[ipnt]*0.5f;
01328
01329 Elo = rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5f;
01330
01331
01332
01333 if( fabs( Elo/exact - 1. ) < 0.001 )
01334 return;
01335
01336 ASSERT( Elo <= exact );
01337
01338 OldEner = rfield.anu[ipnt];
01339
01340
01341 rfield.anu[ipnt] = (realnum)((Ehi + exact)/2.);
01342
01343 rfield.anu[ipnt-1] = (realnum)((exact + Elo)/2.);
01344
01345 rfield.widflx[ipnt] = (realnum)(Ehi - exact);
01346 rfield.widflx[ipnt-1] = (realnum)(exact - Elo);
01347
01348
01349 rfield.anu[ipnt+1] -= (OldEner - rfield.anu[ipnt])/2.f;
01350
01351
01352 ASSERT( rfield.widflx[ipnt-1] > 0. );
01353 ASSERT( rfield.widflx[ipnt] > 0. );
01354 return;
01355 }
01356
01357
01358
01359 STATIC void ipShells(
01360
01361 long int nelem)
01362 {
01363 char chLab[5];
01364 long int
01365 imax,
01366 ion,
01367 nelec,
01368 ns,
01369 nshell;
01370
01371 double thresh=-DBL_MAX;
01372
01373 DEBUG_ENTRY( "ipShells()" );
01374
01375 ASSERT( nelem >= NISO);
01376 ASSERT( nelem < LIMELM );
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387 for( ion=0; ion < nelem; ion++ )
01388 {
01389
01390
01391 strcpy( chLab, elementnames.chElementSym[nelem] );
01392
01393
01394 strcat( chLab, elementnames.chIonStage[ion] );
01395
01396
01397 long int ipISO = nelem-ion;
01398
01399
01400 nelec = ipISO+1;
01401
01402
01403
01404 imax = Heavy.nsShells[nelem][ion];
01405
01406
01407 for( nshell=0; nshell < imax; nshell++ )
01408 {
01409
01410 thresh = (double)(t_ADfA::Inst().ph1(nshell,nelec-1,nelem,0)/EVRYD* 0.9998787);
01411 if( thresh <= 0.1 )
01412 {
01413
01414
01415
01416
01417 opac.ipElement[nelem][ion][nshell][0] = 2;
01418 opac.ipElement[nelem][ion][nshell][1] = 1;
01419 }
01420 else
01421 {
01422
01423
01424
01425
01426 opac.ipElement[nelem][ion][nshell][0] =
01427 ipContEnergy( thresh , chLab );
01428
01429
01430
01431
01432
01433
01434
01435
01436 opac.ipElement[nelem][ion][nshell][1] =
01437 LimitSh(ion+1, nshell+1,nelem+1);
01438 ASSERT( opac.ipElement[nelem][ion][nshell][1] > 0);
01439 }
01440 }
01441
01442 ASSERT( imax > 0 && imax <= 7 );
01443
01444
01445
01446 opac.ipElement[nelem][ion][imax-1][0] =
01447 ipContEnergy(thresh, chLab);
01448
01449
01450 Heavy.ipHeavy[nelem][ion] = opac.ipElement[nelem][ion][imax-1][0];
01451
01452
01453
01454 Heavy.Valence_IP_Ryd[nelem][ion] = thresh;
01455
01456 Heavy.xLyaHeavy[nelem][ion] = 0.;
01457 if( ipISO >= NISO )
01458 {
01459
01460
01461 Heavy.ipLyHeavy[nelem][ion] =
01462 ipLineEnergy(thresh*0.75,chLab , 0);
01463
01464
01465 Heavy.ipBalHeavy[nelem][ion] =
01466 ipLineEnergy(thresh*0.25,chLab , 0);
01467 }
01468 else
01469 {
01470
01471
01472 Heavy.ipLyHeavy[nelem][ion] = -1;
01473 Heavy.ipBalHeavy[nelem][ion] = -1;
01474 }
01475 }
01476
01477
01478
01479 Heavy.nsShells[nelem][nelem] = 1;
01480
01481
01482
01483
01484
01485
01486
01487 opac.ipElement[nelem][nelem][0][0] = ipoint( iso.xIsoLevNIonRyd[ipH_LIKE][nelem][ipH1s] );
01488 ASSERT( opac.ipElement[nelem][nelem][0][0] > 0 );
01489
01490
01491 opac.ipElement[nelem][nelem][0][1] = continuum.KshellLimit;
01492
01493 Heavy.ipHeavy[nelem][nelem] = opac.ipElement[nelem][nelem][0][0];
01494
01495
01496 if( trace.lgTrace && trace.lgPointBug )
01497 {
01498 for( ion=0; ion < (nelem+1); ion++ )
01499 {
01500 fprintf( ioQQQ, "Ion:%3ld%3ld %2.2s%2.2s total shells:%3ld\n",
01501 nelem, ion+1, elementnames.chElementSym[nelem], elementnames.chIonStage[ion]
01502 , Heavy.nsShells[nelem][ion] );
01503 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
01504 {
01505 fprintf( ioQQQ, " shell%3ld %2.2s range eV%10.2e-%8.2e\n",
01506 ns+1, Heavy.chShell[ns], rfield.anu[opac.ipElement[nelem][ion][ns][0]-1]*
01507 EVRYD, rfield.anu[opac.ipElement[nelem][ion][ns][1]-1]*EVRYD );
01508 }
01509 }
01510 }
01511 return;
01512 }
01513
01514
01515 STATIC long LimitSh(long int ion,
01516 long int nshell,
01517 long int nelem)
01518 {
01519 long int LimitSh_v;
01520
01521 DEBUG_ENTRY( "LimitSh()" );
01522
01523
01524
01525
01526 if( nshell == 1 )
01527 {
01528
01529 LimitSh_v = continuum.KshellLimit;
01530
01531 }
01532 else if( nshell == 2 )
01533 {
01534
01535
01536
01537 LimitSh_v = continuum.KshellLimit;
01538
01539 }
01540 else if( nshell == 3 )
01541 {
01542
01543
01544
01545 LimitSh_v = continuum.KshellLimit;
01546
01547 }
01548 else if( nshell == 4 )
01549 {
01550
01551
01552
01553 LimitSh_v = opac.ipElement[nelem-1][ion-1][0][0] - 1;
01554
01555 }
01556 else if( nshell == 5 )
01557 {
01558
01559
01560
01561 LimitSh_v = opac.ipElement[nelem-1][ion-1][0][0] - 1;
01562
01563 }
01564 else if( nshell == 6 )
01565 {
01566
01567
01568
01569 LimitSh_v = opac.ipElement[nelem-1][ion-1][0][0] - 1;
01570
01571 }
01572 else if( nshell == 7 )
01573 {
01574
01575 if( opac.ipElement[nelem-1][ion-1][5][0] < 3 )
01576 {
01577
01578
01579 LimitSh_v = opac.ipElement[nelem-1][ion-1][4][0] -
01580 1;
01581 }
01582 else
01583 {
01584 LimitSh_v = opac.ipElement[nelem-1][ion-1][5][0] -
01585 1;
01586 }
01587
01588 LimitSh_v = opac.ipElement[nelem-1][ion-1][2][0] - 1;
01589
01590 }
01591 else
01592 {
01593 fprintf( ioQQQ, " LimitSh cannot handle nshell as large as%4ld\n",
01594 nshell );
01595 cdEXIT(EXIT_FAILURE);
01596 }
01597 return( LimitSh_v );
01598 }
01599
01600
01601 STATIC void ContBandsCreate(
01602
01603
01604
01605 const char chFile[] )
01606 {
01607 char chLine[FILENAME_PATH_LENGTH_2];
01608 const char* chFilename;
01609 FILE *ioDATA;
01610
01611 bool lgEOL;
01612 long int i,k;
01613
01614
01615
01616 static bool lgCalled=false;
01617
01618 DEBUG_ENTRY( "ContBandsCreate()" );
01619
01620
01621 if( lgCalled )
01622 {
01623
01624 return;
01625 }
01626 lgCalled = true;
01627
01628
01629 chFilename = ( strlen(chFile) == 0 ) ? "bands_continuum.dat" : chFile;
01630
01631
01632 if( trace.lgTrace )
01633 {
01634 fprintf( ioQQQ, " ContBandsCreate opening %s:", chFilename );
01635 }
01636
01637 ioDATA = open_data( chFilename, "r" );
01638
01639
01640 continuum.nContBand = 0;
01641
01642
01643 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01644 {
01645 fprintf( ioQQQ, " ContBandsCreate could not read first line of %s.\n", chFilename );
01646 cdEXIT(EXIT_FAILURE);
01647 }
01648 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
01649 {
01650
01651
01652 if( chLine[0] != '#')
01653 ++continuum.nContBand;
01654 }
01655
01656
01657 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
01658 {
01659 fprintf( ioQQQ, " ContBandsCreate could not rewind %s.\n", chFilename );
01660 cdEXIT(EXIT_FAILURE);
01661 }
01662
01663 continuum.ContBandWavelength = (realnum *)MALLOC(sizeof(realnum)*(unsigned)(continuum.nContBand) );
01664 continuum.chContBandLabels = (char **)MALLOC(sizeof(char *)*(unsigned)(continuum.nContBand) );
01665 continuum.ipContBandLow = (long int *)MALLOC(sizeof(long int)*(unsigned)(continuum.nContBand) );
01666 continuum.ipContBandHi = (long int *)MALLOC(sizeof(long int)*(unsigned)(continuum.nContBand) );
01667
01668
01669 for( i=0; i<continuum.nContBand; ++i )
01670 {
01671
01672 continuum.chContBandLabels[i] = (char *)MALLOC(sizeof(char)*(unsigned)(5) );
01673 }
01674
01675
01676 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01677 {
01678 fprintf( ioQQQ, " ContBandsCreate could not read first line of %s.\n", chFilename );
01679 cdEXIT(EXIT_FAILURE);
01680 }
01681
01682
01683
01684 {
01685 long int m1 , m2 , m3,
01686 myr = 8,
01687 mmo = 12,
01688 mdy = 20;
01689 i = 1;
01690 m1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01691 m2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01692 m3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01693 if( ( m1 != myr ) ||
01694 ( m2 != mmo ) ||
01695 ( m3 != mdy ) )
01696 {
01697 fprintf( ioQQQ,
01698 " ContBandsCreate: the version of the data file %s I found (%li %li %li)is not the current version (%li %li %li).\n",
01699 chFilename ,
01700 m1 , m2 , m3 ,
01701 myr , mmo , mdy );
01702 fprintf( ioQQQ,
01703 " ContBandsCreate: you need to update this file.\n");
01704 cdEXIT(EXIT_FAILURE);
01705 }
01706 }
01707
01708
01709 k = 0;
01710 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
01711 {
01712
01713
01714 if( chLine[0] != '#')
01715 {
01716
01717 strncpy( continuum.chContBandLabels[k] , chLine , 4 );
01718 continuum.chContBandLabels[k][4] = 0;
01719
01720
01721
01722
01723 i = 6;
01724 continuum.ContBandWavelength[k] = (realnum)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01725
01726
01727
01728
01729 continuum.ContBandWavelength[k] *= 1e4f;
01730
01731
01732
01733
01734 double xHi = FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL)*1e4;
01735 double xLow = FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL)*1e4;
01736 if( lgEOL )
01737 {
01738 fprintf( ioQQQ, " There should have been 3 numbers on this band line. Sorry.\n" );
01739 fprintf( ioQQQ, " string==%s==\n" ,chLine );
01740 cdEXIT(EXIT_FAILURE);
01741 }
01742
01743
01744 if( xHi >= xLow )
01745 {
01746 fprintf( ioQQQ, " ContBandWavelength band %li "
01747 "edges are in improper order.\n" ,k);
01748 cdEXIT(EXIT_FAILURE);
01749 }
01750
01751
01752
01753 if( continuum.ContBandWavelength[k] < xHi ||
01754 continuum.ContBandWavelength[k] > xLow )
01755 {
01756 fprintf( ioQQQ, " ContBandWavelength band %li central "
01757 "wavelength not within band.\n" ,k);
01758 cdEXIT(EXIT_FAILURE);
01759 }
01760
01761
01762
01763 continuum.ipContBandHi[k] = ipoint( RYDLAM / xHi );
01764 continuum.ipContBandLow[k] = ipoint( RYDLAM / xLow );
01765
01766
01767
01768
01769
01770
01771 if( trace.lgTrace && trace.lgConBug )
01772 {
01773 if( k==0 )
01774 fprintf( ioQQQ, " ContCreatePointer trace bands\n");
01775 fprintf( ioQQQ,
01776 " band %ld label %s low wl= %.3e low ipnt= %li "
01777 " hi wl= %.3e hi ipnt= %li \n",
01778 k,
01779 continuum.chContBandLabels[k] ,
01780 xLow,
01781 continuum.ipContBandLow[k] ,
01782 xHi,
01783 continuum.ipContBandHi[k] );
01784 }
01785
01786
01787
01788
01789
01790
01791
01792
01793 ++k;
01794 }
01795 }
01796
01797 for( i=0; i<continuum.nContBand; ++i )
01798 {
01799
01800 if( continuum.ContBandWavelength[i] <=0. )
01801 {
01802 fprintf( ioQQQ, " ContBandWavelength band %li has non-positive entry.\n",i );
01803 cdEXIT(EXIT_FAILURE);
01804 }
01805 }
01806
01807 fclose(ioDATA);
01808 return;
01809 }