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