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