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