00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041 #include "cddefines.h"
00042 #include "cddrive.h"
00043 #include "thermal.h"
00044 #include "physconst.h"
00045 #include "doppvel.h"
00046 #include "taulines.h"
00047 #include "dense.h"
00048 #include "rfield.h"
00049 #include "radius.h"
00050 #include "lines_service.h"
00051 #include "ipoint.h"
00052 #include "thirdparty.h"
00053 #include "hydrogenic.h"
00054 #include "lines.h"
00055 #include "rt.h"
00056 #include "trace.h"
00057 #include "save.h"
00058 #include "phycon.h"
00059 #include "atomfeii.h"
00060 #include "iso.h"
00061 #include "pressure.h"
00062 #include "parser.h"
00063
00064
00065 STATIC void FeIIContCreate(double xLamLow , double xLamHigh , long int ncell );
00066
00067
00068
00069
00070
00071 STATIC int FeIIBandsCreate( const char chFile[] );
00072
00073
00074
00075
00076
00077
00078 realnum CS2SMALL = (realnum)1e-5;
00079
00080
00081
00082
00083
00084 STATIC void FeIICollRatesBoltzmann(void);
00085
00086 STATIC void FeIILyaPump(void);
00087
00088
00089
00090
00091
00092
00093
00094 int **ncs1;
00095
00096
00097
00098
00099
00100 #define NPRADDAT 159
00101
00102
00103
00104
00105
00106 realnum **FeII_Bands;
00107
00108
00109
00110
00111
00112
00113 realnum **FeII_Cont;
00114
00115
00116 long int nFeIIBands;
00117
00118
00119 long int nFeIIConBins;
00120
00121
00122
00123
00124 static long int *nnPradDat;
00125
00126
00127
00128
00129 static realnum ***sPradDat;
00130
00131
00132
00133
00134 static double **Fe2SavN;
00135
00136
00137 static double **Fe2A;
00138
00139
00140 static double **Fe2LPump, **Fe2CPump;
00141
00142
00143 realnum *Fe2Energies;
00144
00145
00146 static realnum **Fe2Coll;
00147
00148
00149
00150
00151 static double
00152
00153 *Fe2DepCoef ,
00154
00155 *Fe2LevelPop ,
00156
00157 *Fe2ColDen ,
00158
00159 *FeIIBoltzmann;
00160
00161
00162 static double EnerLyaProf1,
00163 EnerLyaProf4,
00164 PhotOccNumLyaCenter;
00165 static double
00166
00167 *yVector,
00168
00169
00170 **xMatrix ,
00171
00172
00173 *amat;
00174
00175
00176
00177 void FeII_Colden( const char *chLabel )
00178 {
00179 long int n;
00180
00181 DEBUG_ENTRY( "FeII_Colden()" );
00182
00183 if( strcmp(chLabel,"ZERO") == 0 )
00184 {
00185
00186 for( n=0; n < FeII.nFeIILevel_malloc; ++n )
00187 {
00188
00189 Fe2ColDen[n] = 0.;
00190 }
00191 }
00192
00193 else if( strcmp(chLabel,"ADD ") == 0 )
00194 {
00195
00196 for( n=0; n < FeII.nFeIILevel_local; ++n )
00197 {
00198
00199 Fe2ColDen[n] += Fe2LevelPop[n]*radius.drad_x_fillfac;
00200 }
00201 }
00202
00203
00204 else if( strcmp(chLabel,"PRIN") != 0 )
00205 {
00206 fprintf( ioQQQ, " FeII_Colden does not understand the label %s\n",
00207 chLabel );
00208 cdEXIT(EXIT_FAILURE);
00209 }
00210
00211 return;
00212 }
00213
00214
00215
00216
00217
00218
00219 void FeIICreate(void)
00220 {
00221 FILE *ioDATA;
00222
00223 char chLine[FILENAME_PATH_LENGTH_2];
00224
00225 long int i,
00226 ipHi ,
00227 ipLo,
00228 lo,
00229 ihi,
00230 k,
00231 m1,
00232 m2,
00233 m3;
00234
00235 DEBUG_ENTRY( "FeIICreate()" );
00236
00237 if( lgFeIIMalloc )
00238 {
00239
00240
00241 return;
00242 }
00243
00244
00245
00246
00247 lgFeIIMalloc = true;
00248
00249
00250
00251 FeII.nFeIILevel_malloc = FeII.nFeIILevel_local;
00252
00253
00254 Fe2A = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevel_malloc );
00255
00256
00257 for( ipHi=0; ipHi < FeII.nFeIILevel_malloc; ++ipHi )
00258 {
00259 Fe2A[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)FeII.nFeIILevel_malloc);
00260 }
00261
00262
00263 Fe2CPump = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevel_malloc );
00264
00265
00266 Fe2LPump = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevel_malloc );
00267
00268
00269 for( ipHi=0; ipHi < FeII.nFeIILevel_malloc; ++ipHi )
00270 {
00271 Fe2CPump[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)FeII.nFeIILevel_malloc);
00272
00273 Fe2LPump[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)FeII.nFeIILevel_malloc);
00274 }
00275
00276
00277 Fe2Energies = (realnum *)MALLOC(sizeof(realnum)*(unsigned long)FeII.nFeIILevel_malloc );
00278
00279
00280 Fe2Coll = (realnum **)MALLOC(sizeof(realnum *)*(unsigned long)FeII.nFeIILevel_malloc );
00281
00282
00283 for( ipHi=0; ipHi < FeII.nFeIILevel_malloc; ++ipHi )
00284 {
00285 Fe2Coll[ipHi]=(realnum*)MALLOC(sizeof(realnum )*(unsigned long)FeII.nFeIILevel_malloc);
00286 }
00287
00288
00289 xMatrix = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevel_malloc );
00290
00291
00292 for( i=0; i<FeII.nFeIILevel_malloc; ++i )
00293 {
00294 xMatrix[i] = (double *)MALLOC(sizeof(double)*(unsigned long)FeII.nFeIILevel_malloc );
00295 }
00296
00297 amat=(double*)MALLOC( (sizeof(double)*(unsigned long)(FeII.nFeIILevel_malloc*FeII.nFeIILevel_malloc) ));
00298
00299
00300 yVector=(double*)MALLOC( (sizeof(double)*(unsigned long)(FeII.nFeIILevel_malloc) ));
00301
00302
00303 Fe2SavN = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevel_malloc );
00304
00305
00306 for( ipHi=1; ipHi < FeII.nFeIILevel_malloc; ++ipHi )
00307 {
00308 Fe2SavN[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)ipHi);
00309 }
00310
00311
00312 nnPradDat = (long*)MALLOC( (NPRADDAT+1)*sizeof(long) );
00313
00314
00315 Fe2DepCoef = (double*)MALLOC( (unsigned long)FeII.nFeIILevel_malloc*sizeof(double) );
00316
00317
00318 Fe2LevelPop = (double*)MALLOC( (unsigned long)FeII.nFeIILevel_malloc*sizeof(double) );
00319
00320
00321 Fe2ColDen = (double*)MALLOC( (unsigned long)FeII.nFeIILevel_malloc*sizeof(double) );
00322
00323
00324 FeIIBoltzmann = (double*)MALLOC( (unsigned long)FeII.nFeIILevel_malloc*sizeof(double) );
00325
00326
00327
00328
00329 sPradDat = ((realnum ***)MALLOC(NPRADDAT*sizeof(realnum **)));
00330
00331
00332 sPradDat[0] = NULL;
00333 for(ipHi=1; ipHi < NPRADDAT; ipHi++)
00334 {
00335
00336 sPradDat[ipHi] = (realnum **)MALLOC((unsigned long)ipHi*sizeof(realnum *));
00337
00338
00339 for( ipLo=0; ipLo< ipHi; ipLo++ )
00340 {
00341 sPradDat[ipHi][ipLo] = (realnum *)MALLOC(8*sizeof(realnum ));
00342 }
00343 }
00344
00345
00346 for(ipHi=0; ipHi < NPRADDAT; ipHi++)
00347 {
00348 for( ipLo=0; ipLo< ipHi; ipLo++ )
00349 {
00350 for( k=0; k<8; ++k )
00351 {
00352 sPradDat[ipHi][ipLo][k] = -FLT_MAX;
00353 }
00354 }
00355 }
00356
00357
00358 for( ipHi=0; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
00359 {
00360 for( ipLo=0; ipLo < FeII.nFeIILevel_malloc; ipLo++ )
00361 {
00362 FeII.FeIIAul[ipHi][ipLo] = -2.;
00363 for( int k=0; k<8; k++ )
00364 {
00365 FeII.FeIIColl[ipHi][ipLo][k] = -2.;
00366 }
00367 }
00368 FeII.FeIINRGs[ipHi] = -2.;
00369 FeII.FeIISTWT[ipHi] = -2.;
00370 }
00371
00372
00373 ipFe2LevN.reserve(FeII.nFeIILevel_malloc);
00374 ncs1=(int**)MALLOC(sizeof(int *)*(unsigned long)FeII.nFeIILevel_malloc);
00375
00376 for( ipHi=1; ipHi < FeII.nFeIILevel_malloc; ++ipHi )
00377 {
00378 ipFe2LevN.reserve(ipHi,ipHi);
00379 ncs1[ipHi]=(int*)MALLOC(sizeof(int)*(unsigned long)ipHi);
00380 }
00381
00382 ipFe2LevN.alloc();
00383 Fe2LevN.resize(ipFe2LevN.size());
00384 AllTransitions.push_back(Fe2LevN);
00385
00386 int nLine=0;
00387
00388 for( ipHi=1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
00389 {
00390 for( ipLo=0; ipLo < ipHi; ipLo++ )
00391 {
00392 ipFe2LevN[ipHi][ipLo] = nLine;
00393 Fe2LevN[nLine].Junk();
00394 ++nLine;
00395 }
00396 }
00397
00398
00399 Fe2LevN.states()->resize(FeII.nFeIILevel_malloc);
00400
00401 for( ipHi=1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
00402 {
00403
00404 for( ipLo=0; ipLo < ipHi; ipLo++ )
00405 {
00406 TransitionProxy tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
00407 tr.setLo(ipLo);
00408 tr.setHi(ipHi);
00409 tr.AddLine2Stack();
00410 tr.Zero();
00411 }
00412 }
00413
00414 if( trace.lgTrace )
00415 {
00416 fprintf( ioQQQ," FeIICreate opening fe2nn.dat:");
00417 }
00418
00419 ioDATA = open_data( "fe2nn.dat", "r" );
00420
00421 ASSERT( ioDATA !=NULL );
00422
00423
00424
00425
00426
00427
00428 for( i=0; i < 8; i++ )
00429 {
00430 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00431 {
00432 fprintf( ioQQQ, " fe2nn.dat error reading data\n" );
00433 cdEXIT(EXIT_FAILURE);
00434 }
00435
00436
00437 k = 20*i;
00438
00439 ASSERT( k+19 < NPRADDAT+1 );
00440 sscanf( chLine ,
00441 "%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld",
00442 &nnPradDat[k+0], &nnPradDat[k+1], &nnPradDat[k+2], &nnPradDat[k+3], &nnPradDat[k+4],
00443 &nnPradDat[k+5], &nnPradDat[k+6], &nnPradDat[k+7], &nnPradDat[k+8], &nnPradDat[k+9],
00444 &nnPradDat[k+10],&nnPradDat[k+11], &nnPradDat[k+12],&nnPradDat[k+13],&nnPradDat[k+14],
00445 &nnPradDat[k+15],&nnPradDat[k+16], &nnPradDat[k+17],&nnPradDat[k+18],&nnPradDat[k+19]
00446 );
00447 # if !defined(NDEBUG)
00448 for( m1=0; m1<20; ++m1 )
00449 {
00450 ASSERT( nnPradDat[k+m1] >= 0 && nnPradDat[k+m1] <= NFE2LEVN );
00451 }
00452 # endif
00453 }
00454 fclose(ioDATA);
00455
00456
00457 if( trace.lgTrace )
00458 {
00459 fprintf( ioQQQ," FeIICreate opening fe2energies.dat:");
00460 }
00461
00462 ioDATA = open_data( "fe2energies.dat", "r" );
00463
00464
00465 for( ipHi=0; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
00466 {
00467
00468
00469 chLine[0] = '#';
00470 while( chLine[0] == '#' )
00471 {
00472 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00473 {
00474 fprintf( ioQQQ, " fe2energies.dat error reading data\n" );
00475 cdEXIT(EXIT_FAILURE);
00476 }
00477 }
00478
00479
00480 double help;
00481 sscanf( chLine, "%lf", &help );
00482 Fe2Energies[ipHi] = (realnum)help;
00483 FeII.FeIINRGs[ipHi] = Fe2Energies[ipHi];
00484 }
00485 fclose(ioDATA);
00486
00487
00488
00489 if( trace.lgTrace )
00490 fprintf( ioQQQ," FeIICreate opening fe2rad.dat:");
00491
00492 ioDATA = open_data( "fe2rad.dat", "r" );
00493
00494
00495 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00496 {
00497 fprintf( ioQQQ, " fe2rad.dat error reading data\n" );
00498 cdEXIT(EXIT_FAILURE);
00499 }
00500
00501 sscanf( chLine ,"%ld%ld%ld",&lo, &ihi, &m1);
00502 const int nYrFe2Rad=8, nMoFe2Rad=8, nDyFe2Rad=24;
00503 if( lo!=nYrFe2Rad || ihi!=nMoFe2Rad || m1!=nDyFe2Rad )
00504 {
00505 fprintf( ioQQQ, "DISASTER fe2rad.dat has the wrong magic numbers, expected "
00506 "%2i %2i %2i and got %2li %2li %2li\n" ,
00507 nYrFe2Rad, nMoFe2Rad, nDyFe2Rad,
00508 lo, ihi, m1);
00509 cdEXIT(EXIT_FAILURE);
00510 }
00511
00512
00513 for( ipHi=1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
00514 {
00515 for( ipLo=0; ipLo < ipHi; ipLo++ )
00516 {
00517
00518 double save[2];
00519
00520
00521 chLine[0] = '#';
00522 while( chLine[0] == '#' )
00523 {
00524 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00525 {
00526 fprintf( ioQQQ, " fe2nn.dat error reading data\n" );
00527 cdEXIT(EXIT_FAILURE);
00528 }
00529 }
00530
00531
00532 sscanf( chLine ,
00533 "%ld%ld%ld%ld%lf%lf%ld",
00534 &lo, &ihi, &m1, &m2 ,
00535 &save[0], &save[1] , &m3);
00536
00537
00538 const TransitionProxy &tr = Fe2LevN[ipFe2LevN[ihi-1][lo-1]];
00539 (*tr.Lo()).g() = (realnum)m1;
00540 FeII.FeIISTWT[lo-1] = (*tr.Lo()).g();
00541 (*tr.Hi()).g() = (realnum)m2;
00542 FeII.FeIISTWT[ihi-1] = (*tr.Hi()).g();
00543 tr.Emis().Aul() = (realnum)save[0];
00544 FeII.FeIIAul[ihi-1][lo-1] = tr.Emis().Aul();
00545
00546
00547 # define USE_OLD true
00548 if( USE_OLD )
00549 tr.EnergyWN() = (realnum)save[1];
00550 else
00551 tr.EnergyWN() = Fe2Energies[ihi-1]-Fe2Energies[lo-1];
00552
00553
00554 if( fp_equal( tr.Emis().Aul() , (realnum)-1.0f ) )
00555 {
00556
00557 tr.Emis().gf() = 1e-5f;
00558
00559
00560 tr.Emis().Aul() = tr.Emis().gf()*(realnum)TRANS_PROB_CONST*
00561 POW2(tr.EnergyWN())*(*tr.Lo()).g()/(*tr.Hi()).g();
00562 }
00563
00564
00565
00566
00567
00568
00569 ncs1[ihi-1][lo-1] = (int)m3;
00570
00571 }
00572 }
00573 fclose(ioDATA);
00574
00575
00576
00577
00578
00579
00580
00581
00582 if( trace.lgTrace )
00583 {
00584 fprintf( ioQQQ," FeIICreate opening fe2col.dat:");
00585 }
00586
00587 ioDATA = open_data( "fe2col.dat", "r" );
00588
00589 ASSERT( ioDATA !=NULL);
00590 for( ipHi=1; ipHi<NPRADDAT; ++ipHi )
00591 {
00592 for( ipLo=0; ipLo<ipHi; ++ipLo )
00593 {
00594
00595 double save[8];
00596 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00597 {
00598 fprintf( ioQQQ, " fe2col.dat error reading data\n" );
00599 cdEXIT(EXIT_FAILURE);
00600 }
00601 sscanf( chLine ,
00602 "%ld%ld%lf%lf%lf%lf%lf%lf%lf%lf",
00603 &m1, &m2,
00604 &save[0], &save[1] , &save[2],&save[3], &save[4] , &save[5],
00605 &save[6], &save[7]
00606 );
00607 for( k=0; k<8; ++k )
00608 {
00609
00610
00611
00612 sPradDat[m2-1][m1-1][k] = max(CS2SMALL , (realnum)save[k] );
00613
00614 long ipHiFe2 = MAX2( nnPradDat[m2-1] , nnPradDat[m1-1] );
00615 long ipLoFe2 = MIN2( nnPradDat[m2-1] , nnPradDat[m1-1] );
00616 FeII.FeIIColl[ipHiFe2-1][ipLoFe2-1][k] = save[k];
00617 }
00618 }
00619 }
00620
00621 fclose( ioDATA );
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632 for( ipLo=0; ipLo < (FeII.nFeIILevel_malloc - 1); ipLo++ )
00633 {
00634 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
00635 {
00636
00637 const TransitionProxy &tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
00638 ASSERT( tr.EnergyWN() > 0. );
00639 ASSERT( tr.Emis().Aul()> 0. );
00640
00641
00642
00643
00644
00645 tr.WLAng() =
00646 (realnum)(1.0e8/
00647 tr.EnergyWN()/
00648 RefIndex( tr.EnergyWN() ));
00649
00650
00651 tr.Emis().gf() =
00652 (realnum)(tr.Emis().Aul()*(*tr.Hi()).g()/
00653 TRANS_PROB_CONST/POW2(tr.EnergyWN()));
00654
00655 (*tr.Hi()).IonStg() = 2;
00656 (*tr.Hi()).nelem() = 26;
00657
00658
00659
00660
00661
00662
00663
00664
00665 if( ipLo == 0 )
00666 {
00667 tr.Emis().iRedisFun() = FeII.ipRedisFcnResonance;
00668 }
00669 else
00670 {
00671
00672
00673 tr.Emis().iRedisFun() = FeII.ipRedisFcnSubordinate;
00674 }
00675 tr.Emis().phots() = 0.;
00676 tr.Emis().ots() = 0.;
00677 tr.Emis().FracInwd() = 1.;
00678 }
00679 }
00680
00681
00682 if( FeIIBandsCreate("") )
00683 {
00684 fprintf( ioQQQ," failed to open FeII bands file \n");
00685 cdEXIT(EXIT_FAILURE);
00686 }
00687
00688
00689
00690
00691 FeIIContCreate( FeII.feconwlLo , FeII.feconwlHi , FeII.nfe2con );
00692
00693
00694 ASSERT( FeII.nFeIILevel_malloc == FeII.nFeIILevel_local );
00695
00696 return;
00697 }
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717 STATIC bool lgFeIIEverCalled=false;
00718 void FeIILevelPops( void )
00719 {
00720 long int i,
00721 ipHi ,
00722 ipLo ,
00723 n;
00724
00725 bool lgPopNeg;
00726
00727 double
00728 EnergyWN,
00729 pop ,
00730 sum;
00731
00732 int32 info;
00733 int32 ipiv[NFE2LEVN];
00734
00735 DEBUG_ENTRY( "FeIILevelPops()" );
00736
00737 if( trace.lgTrace )
00738 {
00739 fprintf( ioQQQ," FeIILevelPops fe2 pops called\n");
00740 }
00741
00742
00743
00744 if( FeII.lgSimulate )
00745 return;
00746
00747 lgFeIIEverCalled = true;
00748
00749 for( n=0; n<FeII.nFeIILevel_local; ++n)
00750 {
00751 for( ipLo=0; ipLo<FeII.nFeIILevel_local; ++ipLo )
00752 {
00753 Fe2CPump[ipLo][n] = 0.;
00754 Fe2LPump[ipLo][n] = 0.;
00755 Fe2A[ipLo][n] = 0.;
00756 xMatrix[ipLo][n] = 0.;
00757 }
00758 }
00759
00760
00761
00762
00763 FeIICollRatesBoltzmann();
00764
00765
00766 for( n=0; n<FeII.nFeIILevel_local; ++n)
00767 {
00768 for( ipHi=n+1; ipHi<FeII.nFeIILevel_local; ++ipHi )
00769 {
00770 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][n]];
00771
00772 Fe2CPump[n][ipHi] = tr.Emis().pump();
00773
00774
00775 Fe2CPump[ipHi][n] = tr.Emis().pump()*
00776 (*tr.Hi()).g()/(*tr.Lo()).g();
00777
00778
00779 Fe2A[ipHi][n] = tr.Emis().Aul()*(tr.Emis().Pesc() + tr.Emis().Pelec_esc() +
00780 tr.Emis().Pdest());
00781 }
00782 }
00783
00784
00785 FeIILyaPump();
00786
00787
00788
00789
00790
00791
00792
00793
00794 for( n=0; n<FeII.nFeIILevel_local; ++n)
00795 {
00796
00797 for( ipLo=0; ipLo<n; ++ipLo )
00798 {
00799 xMatrix[n][n] = xMatrix[n][n] + Fe2CPump[n][ipLo] + Fe2LPump[n][ipLo]+ Fe2A[n][ipLo] +
00800 Fe2Coll[n][ipLo]*dense.eden;
00801 }
00802
00803 for( ipHi=n+1; ipHi<FeII.nFeIILevel_local; ++ipHi )
00804 {
00805 xMatrix[n][n] = xMatrix[n][n] + Fe2CPump[n][ipHi] + Fe2LPump[n][ipHi] + Fe2Coll[n][ipHi]*dense.eden;
00806 }
00807
00808 for( ipLo=0; ipLo<n; ++ipLo )
00809 {
00810 xMatrix[ipLo][n] = xMatrix[ipLo][n] - Fe2CPump[ipLo][n] - Fe2LPump[ipLo][n] - Fe2Coll[ipLo][n]*dense.eden;
00811 }
00812
00813 for( ipHi=n+1; ipHi<FeII.nFeIILevel_local; ++ipHi )
00814 {
00815 xMatrix[ipHi][n] = xMatrix[ipHi][n] - Fe2CPump[ipHi][n] - Fe2LPump[ipHi][n] - Fe2Coll[ipHi][n]*dense.eden -
00816 Fe2A[ipHi][n];
00817 }
00818
00819
00820 xMatrix[n][0] = 1.0;
00821 }
00822
00823 {
00824
00825 enum {DEBUG_LOC=false};
00826 if( DEBUG_LOC )
00827 {
00828
00829 for( n=0; n<FeII.nFeIILevel_local; ++n)
00830 {
00831
00832
00833 for( ipLo=0; ipLo<FeII.nFeIILevel_local; ++ipLo )
00834 {
00835 fprintf(ioQQQ," %.1e", xMatrix[n][ipLo]);
00836 }
00837 fprintf(ioQQQ,"\n");
00838 }
00839 }
00840 }
00841
00842
00843
00844
00845 yVector[0] = 1.0;
00846 for( n=1; n < FeII.nFeIILevel_local; n++ )
00847 {
00848 yVector[n] = 0.0;
00849 }
00850
00851
00852
00853 # ifdef AMAT
00854 # undef AMAT
00855 # endif
00856 # define AMAT(I_,J_) (*(amat+(I_)*FeII.nFeIILevel_local+(J_)))
00857
00858
00859 for( ipHi=0; ipHi < FeII.nFeIILevel_local; ipHi++ )
00860 {
00861 for( i=0; i < FeII.nFeIILevel_local; i++ )
00862 {
00863 AMAT(i,ipHi) = xMatrix[i][ipHi];
00864 }
00865 }
00866
00867 info = 0;
00868
00869
00870 getrf_wrapper(FeII.nFeIILevel_local, FeII.nFeIILevel_local, amat, FeII.nFeIILevel_local, ipiv, &info);
00871 getrs_wrapper('N', FeII.nFeIILevel_local, 1, amat, FeII.nFeIILevel_local, ipiv, yVector, FeII.nFeIILevel_local, &info);
00872
00873 if( info != 0 )
00874 {
00875 fprintf( ioQQQ, "DISASTER FeIILevelPops: dgetrs finds singular or ill-conditioned matrix\n" );
00876 cdEXIT(EXIT_FAILURE);
00877 }
00878
00879
00880
00881
00882 lgPopNeg = false;
00883
00884 for( ipLo=0; ipLo < FeII.nFeIILevel_local; ipLo++ )
00885 {
00886 if(yVector[ipLo] < 0. )
00887 {
00888 lgPopNeg = true;
00889 fprintf(ioQQQ,"PROBLEM FeIILevelPops finds non-positive level population, level is %ld pop is %g\n",
00890 ipLo , yVector[ipLo] );
00891 }
00892
00893 Fe2LevelPop[ipLo] = yVector[ipLo] * dense.xIonDense[ipIRON][1];
00894 }
00895 if( lgPopNeg )
00896 {
00897
00898
00899 fprintf(ioQQQ , "PROBLEM FeIILevelPops exits with negative level populations.\n");
00900 }
00901
00902
00903
00904
00905
00906 for( ipLo=FeII.nFeIILevel_local; ipLo < FeII.nFeIILevel_malloc; ++ipLo )
00907 {
00908 Fe2LevelPop[ipLo] = 0.;
00909 }
00910
00911
00912
00913
00914 for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
00915 {
00916
00917
00918 for( ipHi=ipLo+1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
00919 {
00920
00921
00922
00923
00924 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
00925 tr.Emis().PopOpc() = (Fe2LevelPop[ipLo] -
00926 Fe2LevelPop[ipHi]*(*tr.Lo()).g()/(*tr.Hi()).g());
00927
00928 (*tr.Lo()).Pop() = Fe2LevelPop[ipLo];
00929
00930 (*tr.Hi()).Pop() = Fe2LevelPop[ipHi];
00931
00932 tr.Emis().phots() = Fe2LevelPop[ipHi]*
00933 tr.Emis().Aul()*(tr.Emis().Pesc() + tr.Emis().Pelec_esc() );
00934
00935 tr.Emis().xIntensity() = tr.Emis().phots()*
00936 tr.EnergyErg();
00937
00938
00939
00940 tr.Emis().ColOvTot() = (realnum)(Fe2Coll[ipLo][ipHi]*dense.eden /
00941 SDIV( Fe2Coll[ipLo][ipHi]*dense.eden + Fe2CPump[ipLo][ipHi] + Fe2LPump[ipLo][ipHi] ) );
00942 }
00943 }
00944
00945
00946
00947 if( FeII.lgFeIILargeOn )
00948 {
00949
00950 hydro.dstfe2lya = 0.;
00951 EnergyWN = 0.;
00952
00953 for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
00954 {
00955 for( ipHi=ipLo+1; ipHi < FeII.nFeIILevel_local; ipHi++ )
00956 {
00957 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
00958 EnergyWN += Fe2LPump[ipLo][ipHi] + Fe2LPump[ipHi][ipLo];
00959 hydro.dstfe2lya += (realnum)(
00960 (*tr.Lo()).Pop()*Fe2LPump[ipLo][ipHi] -
00961 (*tr.Hi()).Pop()*Fe2LPump[ipHi][ipLo]);
00962 }
00963 }
00964
00965
00966 pop = iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop();
00967 if( pop > SMALLFLOAT )
00968 {
00969 hydro.dstfe2lya /= (realnum)(pop * iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Aul());
00970 }
00971 else
00972 {
00973 hydro.dstfe2lya = 0.;
00974 }
00975
00976
00977 }
00978
00979 {
00980 enum {DEBUG_LOC=false};
00981 if( DEBUG_LOC)
00982 {
00983
00984 fprintf(ioQQQ," sum all %.1e dest rate%.1e escR= %.1e\n",
00985 EnergyWN,hydro.dstfe2lya,
00986 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pesc());
00987
00988 }
00989 }
00990
00991
00992
00993
00994 Fe2DepCoef[0] = 1.0;
00995 sum = 1.0;
00996 for( i=1; i < FeII.nFeIILevel_local; i++ )
00997 {
00998
00999 Fe2DepCoef[i] = Fe2DepCoef[0]*FeIIBoltzmann[i]*
01000 (*Fe2LevN[ipFe2LevN[i][0]].Hi()).g()/ (*Fe2LevN[ipFe2LevN[1][0]].Lo()).g();
01001
01002 sum += Fe2DepCoef[i];
01003 }
01004
01005
01006 for( i=0; i < FeII.nFeIILevel_local; i++ )
01007 {
01008
01009
01010 Fe2DepCoef[i] /= sum;
01011
01012
01013 if( Fe2DepCoef[i]>SMALLFLOAT )
01014 {
01015 Fe2DepCoef[i] = yVector[i]/Fe2DepCoef[i];
01016 }
01017 else
01018 {
01019 Fe2DepCoef[i] = 0.;
01020 }
01021 }
01022
01023
01024
01025 FeII.Fe2_large_cool = 0.0f;
01026
01027 FeII.Fe2_large_heat = 0.f;
01028
01029 FeII.ddT_Fe2_large_cool = 0.0f;
01030
01031
01032 for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
01033 {
01034 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_local; ipHi++ )
01035 {
01036 double OneLine;
01037 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
01038
01039
01040 OneLine = (Fe2Coll[ipLo][ipHi]*Fe2LevelPop[ipLo] - Fe2Coll[ipHi][ipLo]*Fe2LevelPop[ipHi])*
01041 dense.eden*tr.EnergyErg();
01042
01043
01044 FeII.Fe2_large_cool += (realnum)MAX2(0., OneLine);
01045
01046
01047 FeII.Fe2_large_heat += (realnum)MAX2(0., -OneLine);
01048
01049
01050 if( OneLine >= 0. )
01051 {
01052
01053 FeII.ddT_Fe2_large_cool += (realnum)OneLine*
01054 (Fe2LevN[ipFe2LevN[ipHi][0]].EnergyK()*thermal.tsq1 - thermal.halfte);
01055 }
01056 else
01057 {
01058
01059 FeII.ddT_Fe2_large_cool -= (realnum)OneLine*thermal.halfte;
01060 }
01061 }
01062 }
01063
01064 return;
01065 }
01066
01067
01068
01069
01070
01071
01072 STATIC void FeIICollRatesBoltzmann(void)
01073 {
01074
01075 static double OldTemp = -1.;
01076 long int i,
01077 ipLo ,
01078 ipHi,
01079 ipTrim;
01080 realnum ag,
01081 cg,
01082 dg,
01083 gb,
01084 y;
01085 realnum FracLowTe , FracHighTe;
01086 static realnum tt[8]={1e3f,3e3f,5e3f,7e3f,1e4f,12e3f,15e3f,2e4f};
01087 long int ipTemp,
01088 ipTempp1 ,
01089 ipLoFe2,
01090 ipHiFe2;
01091 static long int nFeII_old = -1;
01092
01093 DEBUG_ENTRY( "FeIICollRatesBoltzmann()" );
01094
01095
01096
01097 if( fp_equal( phycon.te, OldTemp ) && FeII.nFeIILevel_local == nFeII_old )
01098 {
01099
01100 return;
01101 }
01102 OldTemp = phycon.te;
01103 nFeII_old = FeII.nFeIILevel_local;
01104
01105
01106
01107
01108 for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
01109 {
01110 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_local; ipHi++ )
01111 {
01112
01113 if( ncs1[ipHi][ipLo] == 3 )
01114 {
01115
01116 ag = 0.15f;
01117 cg = 0.f;
01118 dg = 0.f;
01119 }
01120
01121 else if( ncs1[ipHi][ipLo] == 1 )
01122 {
01123 ag = 0.6f;
01124 cg = 0.f;
01125 dg = 0.28f;
01126 }
01127
01128 else if( ncs1[ipHi][ipLo] == 2 )
01129 {
01130 ag = 0.0f;
01131 cg = 0.1f;
01132 dg = 0.0f;
01133 }
01134 else
01135 {
01136
01137 ag = 0.0f;
01138 cg = 0.1f;
01139 dg = 0.0f;
01140 fprintf(ioQQQ,">>>PROBLEM impossible ncs1 in FeIILevelPops\n");
01141 }
01142
01143
01144 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
01145 y = tr.EnergyWN()/ (realnum)phycon.te_wn;
01146
01147 gb = (realnum)(ag + (-cg*POW2(y) + dg)*(log(1.0+1.0/y) - 0.4/
01148 POW2(y + 1.0)) + cg*y);
01149
01150 tr.Coll().col_str() = 23.861f*1e5f*gb*
01151 tr.Emis().Aul()*(*tr.Hi()).g()/
01152 POW3(tr.EnergyWN());
01153
01154
01155 ASSERT( tr.Coll().col_str() > 0.);
01156
01157
01158
01159 tr.Coll().col_str() = MAX2( CS2SMALL, tr.Coll().col_str());
01160
01161
01162
01163
01164
01165
01166
01167 }
01168 }
01169
01170
01171
01172
01173 if( phycon.te <= tt[0] )
01174 {
01175
01176
01177
01178
01179 ipTemp = 0;
01180 ipTempp1 = 0;
01181 FracHighTe = 0.;
01182 }
01183 else if( phycon.te > tt[7] )
01184 {
01185
01186
01187
01188
01189 ipTemp = 7;
01190 ipTempp1 = 7;
01191 FracHighTe = 0.;
01192 }
01193 else
01194 {
01195
01196 ipTemp = -1;
01197 for( i=0; i < 8-1; i++ )
01198 {
01199 if( phycon.te <= tt[i+1] )
01200 {
01201 ipTemp = i;
01202 break;
01203 }
01204
01205 }
01206
01207 if( ipTemp < 0 )
01208 {
01209 fprintf( ioQQQ, " Insanity while looking for temperature in coll str array, te=%g.\n",
01210 phycon.te );
01211 cdEXIT(EXIT_FAILURE);
01212 }
01213
01214
01215 ipTemp = i;
01216 ipTempp1 = i+1;
01217 FracHighTe = ((realnum)phycon.te - tt[ipTemp])/(tt[ipTempp1] - tt[ipTemp]);
01218 }
01219
01220
01221
01222 FracLowTe = 1.f - FracHighTe;
01223
01224 for( ipHi=1; ipHi < NPRADDAT; ipHi++ )
01225 {
01226 for( ipLo=0; ipLo < ipHi; ipLo++ )
01227 {
01228
01229
01230 ipHiFe2 = MAX2( nnPradDat[ipHi] , nnPradDat[ipLo] );
01231 ipLoFe2 = MIN2( nnPradDat[ipHi] , nnPradDat[ipLo] );
01232 ASSERT( ipHiFe2-1 < NFE2LEVN );
01233 ASSERT( ipHiFe2-1 >= 0 );
01234 ASSERT( ipLoFe2-1 < NFE2LEVN );
01235 ASSERT( ipLoFe2-1 >= 0 );
01236
01237
01238 if( ipHiFe2-1 < FeII.nFeIILevel_local )
01239 {
01240 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHiFe2-1][ipLoFe2-1]];
01241
01242
01243 tr.Coll().col_str() =
01244 sPradDat[ipHi][ipLo][ipTemp] * FracLowTe +
01245 sPradDat[ipHi][ipLo][ipTempp1] * FracHighTe;
01246
01247
01248
01249 ASSERT( tr.Coll().col_str() > 0. );
01250 }
01251 }
01252 }
01253
01254
01255 FeIIBoltzmann[0] = 1.0;
01256 for( ipHi=1; ipHi < FeII.nFeIILevel_local; ipHi++ )
01257 {
01258
01259
01260
01261
01262 FeIIBoltzmann[ipHi] = dsexp( Fe2LevN[ipFe2LevN[ipHi][0]].EnergyWN()/phycon.te_wn );
01263 }
01264
01265
01266 ipTrim = FeII.nFeIILevel_local - 1;
01267 while( FeIIBoltzmann[ipTrim] == 0. && ipTrim > 0 )
01268 {
01269 --ipTrim;
01270 }
01271
01272
01273
01274
01275 ASSERT( ipTrim > 0 );
01276
01277
01278
01279 if( ipTrim+1 < FeII.nFeIILevel_local )
01280 {
01281
01282 for( ipLo=0; ipLo<FeII.nFeIILevel_local; ++ipLo )
01283 {
01284 for( ipHi=ipTrim; ipHi<FeII.nFeIILevel_local; ++ipHi )
01285 {
01286 Fe2Coll[ipLo][ipHi] = 0.;
01287 Fe2Coll[ipHi][ipLo] = 0.;
01288 }
01289 }
01290 }
01291 FeII.nFeIILevel_local = ipTrim+1;
01292
01293
01294
01295 {
01296 enum {DEBUG_LOC=false};
01297 if( DEBUG_LOC)
01298 {
01299 for( ipLo=0; ipLo < 52; ipLo++ )
01300 {
01301 fprintf(ioQQQ,"%e %e\n",
01302 Fe2LevN[ipFe2LevN[51][ipLo]].Coll().col_str(),Fe2LevN[ipFe2LevN[52][ipLo]].Coll().col_str());
01303 }
01304 cdEXIT(EXIT_FAILURE);
01305 }
01306 }
01307
01308
01309 for( ipLo=0; ipLo<FeII.nFeIILevel_local; ++ipLo)
01310 {
01311 for( ipHi=ipLo+1; ipHi<FeII.nFeIILevel_local; ++ipHi )
01312 {
01313 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
01314
01315
01316
01317 Fe2Coll[ipHi][ipLo] = (realnum)(COLL_CONST/phycon.sqrte*tr.Coll().col_str()/
01318 (*tr.Hi()).g());
01319
01320
01321
01322
01323
01324 Fe2Coll[ipLo][ipHi] = (realnum)(Fe2Coll[ipHi][ipLo]*FeIIBoltzmann[ipHi]/FeIIBoltzmann[ipLo]*
01325 (*tr.Hi()).g()/(*tr.Lo()).g());
01326
01327 }
01328 }
01329
01330 return;
01331 }
01332
01333
01334
01335
01336
01337
01338
01339
01340 void FeIIPrint(void)
01341 {
01342
01343 DEBUG_ENTRY( "FeIIPrint()" );
01344
01345 return;
01346 }
01347
01348
01349
01350
01351
01352
01353
01354
01355 double FeIISumBand(
01356
01357 realnum wl1,
01358 realnum wl2,
01359 double *SumBandInward )
01360 {
01361 long int ipHi,
01362 ipLo;
01363 double SumBandFe2_v;
01364
01365 DEBUG_ENTRY( "FeIISumBand()" );
01366
01367
01368 SumBandFe2_v = 0.;
01369 *SumBandInward = 0.;
01370 if( dense.xIonDense[ipIRON][1] > SMALLFLOAT )
01371 {
01372
01373 ASSERT( wl2 > wl1 );
01374 for( ipHi=1; ipHi < FeII.nFeIILevel_local; ++ipHi )
01375 {
01376 for( ipLo=0; ipLo < ipHi; ++ipLo )
01377 {
01378 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
01379 if( tr.WLAng() >= wl1 &&
01380 tr.WLAng() < wl2 )
01381 {
01382 SumBandFe2_v += tr.Emis().xIntensity();
01383 *SumBandInward += tr.Emis().xIntensity() *
01384 tr.Emis().FracInwd();
01385 }
01386 }
01387 }
01388 }
01389 return( SumBandFe2_v );
01390 }
01391
01392
01393
01394
01395
01396 void FeII_RT_TauInc(void)
01397 {
01398 long int ipHi,
01399 ipLo;
01400
01401 DEBUG_ENTRY( "FeII_RT_TauInc()" );
01402
01403 for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
01404 {
01405
01406
01407
01408
01409 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
01410 {
01411
01412 const TransitionProxy &tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
01413 if( tr.ipCont() > 0 )
01414 RT_line_one_tauinc( tr, -8, -8, ipHi, ipLo, GetDopplerWidth(dense.AtomicWeight[ipIRON]) );
01415 }
01416 }
01417
01418 return;
01419 }
01420
01421
01422
01423
01424
01425 void FeII_RT_tau_reset(void)
01426 {
01427 long int ipHi,
01428 ipLo;
01429
01430 DEBUG_ENTRY( "FeII_RT_tau_reset()" );
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441 for( ipHi=1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
01442 {
01443 for( ipLo=0; ipLo < ipHi; ipLo++ )
01444 {
01445 RT_line_one_tau_reset( Fe2LevN[ipFe2LevN[ipHi][ipLo]] );
01446 }
01447 }
01448
01449 return;
01450 }
01451
01452
01453
01454
01455
01456 void FeIIPoint(void)
01457 {
01458 long int ipHi,
01459 ip,
01460 ipLo;
01461
01462 DEBUG_ENTRY( "FeIIPoint()" );
01463
01464
01465
01466 for( ipLo=0; ipLo < FeII.nFeIILevel_malloc-1; ipLo++ )
01467 {
01468 for( ipHi=ipLo+1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
01469 {
01470
01471
01472 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
01473 if( fabs(tr.Emis().Aul()- 1e-5 ) > 1e-8 )
01474 {
01475 ip = ipoint(tr.EnergyRyd());
01476 tr.ipCont() = ip;
01477
01478
01479 if( strcmp( rfield.chLineLabel[ip-1], " " ) == 0 )
01480 strcpy( rfield.chLineLabel[ip-1], "FeII" );
01481
01482 tr.Emis().ipFine() = ipFineCont( tr.EnergyRyd());
01483 }
01484 else
01485 {
01486 tr.ipCont() = -1;
01487 tr.Emis().ipFine() = -1;
01488 }
01489
01490 tr.Emis().dampXvel() =
01491 (realnum)(tr.Emis().Aul()/
01492 tr.EnergyWN()/PI4);
01493
01494
01495 tr.Emis().opacity() =
01496 (realnum)abscf(tr.Emis().gf(),
01497 tr.EnergyWN(),
01498 (*tr.Lo()).g());
01499 }
01500 }
01501
01502 return;
01503 }
01504
01505
01506
01507
01508
01509 void FeIIAccel(double *fe2drive)
01510 {
01511 long int ipHi,
01512 ipLo;
01513
01514 DEBUG_ENTRY( "FeIIAccel()" );
01515
01516
01517
01518
01519
01520 *fe2drive = 0.;
01521 for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
01522 {
01523 for( ipHi=ipLo+1; ipHi < FeII.nFeIILevel_local; ipHi++ )
01524 {
01525 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
01526 *fe2drive += tr.Emis().pump()*
01527 tr.EnergyErg()*tr.Emis().PopOpc();
01528 }
01529 }
01530
01531 return;
01532 }
01533
01534
01535
01536
01537
01538 void FeII_RT_Make( void )
01539 {
01540 long int ipHi,
01541 ipLo;
01542
01543 DEBUG_ENTRY( "FeII_RT_Make()" );
01544
01545 if( trace.lgTrace )
01546 {
01547 fprintf( ioQQQ," FeII_RT_Make called\n");
01548 }
01549
01550
01551 for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
01552 {
01553
01554 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
01555 {
01556
01557 const TransitionProxy &tr=Fe2LevN[ipFe2LevN[ipHi][ipLo]];
01558 if( tr.ipCont() > 0 )
01559 {
01560
01561
01562 RT_line_one( tr, true, 0.f, GetDopplerWidth(dense.AtomicWeight[ipIRON]) );
01563 }
01564 }
01565 }
01566
01567 return;
01568 }
01569
01570
01571
01572
01573
01574 void FeIIAddLines( void )
01575 {
01576
01577 DEBUG_ENTRY( "FeIIAddLines()" );
01578
01579
01580
01581
01582
01583
01584
01585
01586
01587 if( LineSave.ipass == 0 )
01588 {
01589
01590 for( long ipLo=0; ipLo < (FeII.nFeIILevel_malloc - 1); ipLo++ )
01591 {
01592 for( long ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
01593 {
01594 Fe2SavN[ipHi][ipLo] = 0.;
01595 }
01596 }
01597 }
01598
01599
01600 else if( LineSave.ipass == 1 )
01601 {
01602
01603 for( long ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
01604 {
01605 for( long ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_local; ipHi++ )
01606 {
01607 Fe2SavN[ipHi][ipLo] +=
01608 radius.dVeffAper*Fe2LevN[ipFe2LevN[ipHi][ipLo]].Emis().xIntensity();
01609 }
01610 }
01611 }
01612
01613 {
01614 enum {DEBUG_LOC=false};
01615 if( DEBUG_LOC )
01616 {
01617 fprintf(ioQQQ," 69-1\t%li\t%e\n", nzone , Fe2SavN[68][0] );
01618 }
01619 }
01620
01621 return;
01622 }
01623
01624
01625 void FeIIPunchLevels(
01626
01627 FILE *ioPUN )
01628 {
01629
01630 long int ipHi;
01631
01632 DEBUG_ENTRY( "FeIIPunchLevels()" );
01633
01634 fprintf(ioPUN , "%.2f\t%li\n", 0., (long)(*Fe2LevN[ipFe2LevN[1][0]].Lo()).g() );
01635 for( ipHi=1; ipHi < FeII.nFeIILevel_malloc; ++ipHi )
01636 {
01637 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][0]];
01638 fprintf(ioPUN , "%.2f\t%li\n",
01639 tr.EnergyWN() ,
01640 (long)(*tr.Hi()).g());
01641 }
01642
01643 return;
01644 }
01645
01646
01647 void FeIIPunchColden(
01648
01649 FILE *ioPUN )
01650 {
01651
01652 long int n;
01653
01654 DEBUG_ENTRY( "FeIIPunchColden()" );
01655
01656 fprintf(ioPUN , "%.2f\t%li\t%.3e\n", 0., (long)(*Fe2LevN[ipFe2LevN[1][0]].Lo()).g() , Fe2ColDen[0]);
01657 for( n=1; n < FeII.nFeIILevel_malloc; ++n )
01658 {
01659 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[n][0]];
01660 fprintf(ioPUN , "%.2f\t%li\t%.3e\n",
01661 tr.EnergyWN() ,
01662 (long)(*tr.Hi()).g(),
01663 Fe2ColDen[n]);
01664 }
01665
01666 return;
01667 }
01668
01669
01670
01671
01672
01673 void FeIIPunchOpticalDepth(
01674
01675 FILE *ioPUN )
01676 {
01677 long int
01678 ipHi,
01679 ipLo;
01680
01681 DEBUG_ENTRY( "FeIIPunchOpticalDepth()" );
01682
01683
01684
01685
01686 for( ipLo=0; ipLo < (FeII.nFeIILevel_malloc - 1); ipLo++ )
01687 {
01688 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
01689 {
01690
01691
01692
01693
01694 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
01695 fprintf( ioPUN, "%ld\t%ld\t%.2f\t%.2e\n",
01696 ipHi+1,
01697 ipLo+1,
01698 tr.WLAng() ,
01699 tr.Emis().TauIn() );
01700 }
01701 }
01702
01703 return;
01704 }
01705
01706
01707
01708
01709 void FeIISaveLines(
01710
01711 FILE *ioPUN )
01712 {
01713 long int MaseHi,
01714 MaseLow,
01715 ipHi,
01716 ipLo;
01717 double hbeta, absint , renorm;
01718
01719 realnum TauMase,
01720 thresh;
01721
01722 DEBUG_ENTRY( "FeIISaveLines()" );
01723
01724
01725
01726
01727
01728 if( LineSv[LineSave.ipNormWavL].SumLine[0] > 0. )
01729 renorm = LineSave.ScaleNormLine/LineSv[LineSave.ipNormWavL].SumLine[0];
01730 else
01731 renorm = 1.;
01732
01733 fprintf( ioPUN, " up low log I, I/I(LineSave), Tau\n" );
01734
01735
01736 MaseLow = -1;
01737 MaseHi = -1;
01738 TauMase = 0.f;
01739 for( ipLo=0; ipLo < (FeII.nFeIILevel_malloc - 1); ipLo++ )
01740 {
01741 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
01742 {
01743 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
01744 if( tr.Emis().TauIn() < TauMase )
01745 {
01746 TauMase = tr.Emis().TauIn();
01747 MaseLow = ipLo;
01748 MaseHi = ipHi;
01749 }
01750 }
01751 }
01752
01753 if( TauMase < 0.f )
01754 {
01755 fprintf( ioPUN, " Most negative optical depth was %4ld%4ld%10.2e\n",
01756 MaseHi, MaseLow, TauMase );
01757 }
01758
01759
01760 if( cdLine("TOTL", 4861.36f , &hbeta , &absint)<=0 )
01761 {
01762 fprintf( ioQQQ, " FeIILevelPops could not find Hbeta with cdLine.\n" );
01763 cdEXIT(EXIT_FAILURE);
01764 }
01765
01766 fprintf( ioPUN, "Hbeta=%7.3f %10.2e\n",
01767 absint ,
01768 hbeta );
01769
01770 if( renorm > SMALLFLOAT )
01771
01772
01773 thresh = FeII.fe2thresh/(realnum)renorm;
01774 else
01775 thresh = 0.f;
01776
01777
01778
01779 for( ipLo=0; ipLo < (FeII.nFeIILevel_malloc - 1); ipLo++ )
01780 {
01781 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
01782 {
01783
01784
01785
01786
01787 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
01788 if( (Fe2SavN[ipHi][ipLo] > thresh &&
01789 tr.EnergyWN() > FeII.fe2ener[0]) &&
01790 tr.EnergyWN() < FeII.fe2ener[1] )
01791 {
01792 if( FeII.lgShortFe2 )
01793 {
01794
01795
01796 fprintf( ioPUN, "%ld\t%ld\t%.2f\t%.3f\n",
01797 ipHi+1,
01798 ipLo+1,
01799 tr.WLAng() ,
01800 log10(MAX2(1e-37,Fe2SavN[ipHi][ipLo])) + radius.Conv2PrtInten );
01801 }
01802 else
01803 {
01804
01805 fprintf( ioPUN, "%ld\t%ld\t%.2f\t%.3f\t%.2e\t%.2e\n",
01806 ipHi+1,
01807 ipLo+1,
01808 tr.WLAng() ,
01809 log10(MAX2(1e-37,Fe2SavN[ipHi][ipLo])) + radius.Conv2PrtInten,
01810 Fe2SavN[ipHi][ipLo]*renorm ,
01811 tr.Emis().TauIn() );
01812 }
01813 }
01814 }
01815 }
01816
01817 return;
01818 }
01819
01820
01821
01822
01823
01824
01825 void FeII_LineZero(void)
01826 {
01827 long int ipHi,
01828 ipLo;
01829
01830 DEBUG_ENTRY( "FeII_LineZero()" );
01831
01832
01833
01834
01835
01836 for( ipLo=0; ipLo < (FeII.nFeIILevel_malloc - 1); ipLo++ )
01837 {
01838 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
01839 {
01840
01841 Fe2LevN[ipFe2LevN[ipHi][ipLo]].Zero();
01842 }
01843 }
01844
01845 return;
01846 }
01847
01848
01849
01850
01851 void FeIIIntenZero(void)
01852 {
01853 long int ipHi,
01854 ipLo;
01855
01856 DEBUG_ENTRY( "FeIIIntenZero()" );
01857
01858
01859
01860
01861 Fe2LevelPop[0] = 0.;
01862 for( ipHi=1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
01863 {
01864
01865 Fe2LevelPop[ipHi] = 0.;
01866 for( ipLo=0; ipLo < ipHi; ipLo++ )
01867 {
01868 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
01869
01870
01871 tr.Emis().PopOpc() = 0.;
01872
01873
01874 (*tr.Lo()).Pop() = 0.;
01875
01876
01877 (*tr.Hi()).Pop() = 0.;
01878
01879
01880 tr.Coll().cool() = 0.;
01881 tr.Coll().heat() = 0.;
01882
01883
01884 tr.Emis().xIntensity() = 0.;
01885
01886 tr.Emis().phots() = 0.;
01887 tr.Emis().ots() = 0.;
01888 tr.Emis().ColOvTot() = 0.;
01889 }
01890 }
01891
01892 (*TauLines[ipTuv3].Lo()).Pop() = 0;
01893 (*TauLines[ipTuv3].Hi()).Pop() = 0;
01894 TauLines[ipTuv3].Emis().PopOpc() = 0;
01895 TauLines[ipTuv3].Emis().phots() = 0;
01896 TauLines[ipTuv3].Emis().xIntensity() = 0;
01897
01898 (*TauLines[ipTr48].Lo()).Pop() = 0;
01899 (*TauLines[ipTr48].Hi()).Pop() = 0;
01900 TauLines[ipTr48].Emis().PopOpc() = 0;
01901 TauLines[ipTr48].Emis().phots() = 0;
01902 TauLines[ipTr48].Emis().xIntensity() = 0;
01903
01904 FeII.for7 = 0;
01905
01906 (*TauLines[ipTFe16].Lo()).Pop() = 0;
01907 (*TauLines[ipTFe16].Hi()).Pop() = 0;
01908 TauLines[ipTFe16].Emis().PopOpc() = 0;
01909 TauLines[ipTFe16].Emis().phots() = 0;
01910 TauLines[ipTFe16].Emis().xIntensity() = 0;
01911
01912 (*TauLines[ipTFe26].Lo()).Pop() = 0;
01913 (*TauLines[ipTFe26].Hi()).Pop() = 0;
01914 TauLines[ipTFe26].Emis().PopOpc() = 0;
01915 TauLines[ipTFe26].Emis().phots() = 0;
01916 TauLines[ipTFe26].Emis().xIntensity() = 0;
01917
01918 (*TauLines[ipTFe34].Lo()).Pop() = 0;
01919 (*TauLines[ipTFe34].Hi()).Pop() = 0;
01920 TauLines[ipTFe34].Emis().PopOpc() = 0;
01921 TauLines[ipTFe34].Emis().phots() = 0;
01922 TauLines[ipTFe34].Emis().xIntensity() = 0;
01923
01924 (*TauLines[ipTFe35].Lo()).Pop() = 0;
01925 (*TauLines[ipTFe35].Hi()).Pop() = 0;
01926 TauLines[ipTFe35].Emis().PopOpc() = 0;
01927 TauLines[ipTFe35].Emis().phots() = 0;
01928 TauLines[ipTFe35].Emis().xIntensity() = 0;
01929
01930 (*TauLines[ipTFe46].Lo()).Pop() = 0;
01931 (*TauLines[ipTFe46].Hi()).Pop() = 0;
01932 TauLines[ipTFe46].Emis().PopOpc() = 0;
01933 TauLines[ipTFe46].Emis().phots() = 0;
01934 TauLines[ipTFe46].Emis().xIntensity() = 0;
01935
01936 (*TauLines[ipTFe56].Lo()).Pop() = 0;
01937 (*TauLines[ipTFe56].Hi()).Pop() = 0;
01938 TauLines[ipTFe56].Emis().PopOpc() = 0;
01939 TauLines[ipTFe56].Emis().phots() = 0;
01940 TauLines[ipTFe56].Emis().xIntensity() = 0;
01941
01942 (*TauLines[ipT191].Lo()).Pop() = 0;
01943 (*TauLines[ipT191].Hi()).Pop() = 0;
01944 TauLines[ipT191].Emis().PopOpc() = 0;
01945 TauLines[ipT191].Emis().phots() = 0;
01946 TauLines[ipT191].Emis().xIntensity() = 0;
01947
01948 return;
01949 }
01950
01951
01952
01953
01954 void FeIIPunData(
01955
01956 FILE* ioPUN ,
01957
01958 bool lgDoAll )
01959 {
01960 long int ipLo , ipHi;
01961
01962 DEBUG_ENTRY( "FeIIPunData()" );
01963
01964 if( lgDoAll )
01965 {
01966 fprintf( ioQQQ,
01967 " FeIIPunData ALL option not implemented yet 1\n" );
01968 cdEXIT(EXIT_FAILURE);
01969 }
01970 else if( lgFeIIEverCalled )
01971 {
01972 long int nSkip=0, limit=MIN2(64, FeII.nFeIILevel_local);
01973
01974
01975
01976
01977 bool lgPrint = true;
01978 for( ipHi=1; ipHi<limit; ++ipHi )
01979 {
01980 for( ipLo=0; ipLo<ipHi; ++ipLo )
01981 {
01982
01983 Save1LineData( Fe2LevN[ipFe2LevN[ipHi][ipLo]] , ioPUN, false , lgPrint);
01984 }
01985 }
01986 fprintf( ioPUN , "\n");
01987
01988 if( limit==64 )
01989 {
01990
01991
01992
01993 for( ipHi=limit; ipHi<FeII.nFeIILevel_local; ++ipHi )
01994 {
01995
01996
01997 for( ipLo=0; ipLo<ipHi; ++ipLo )
01998 {
01999
02000
02001
02002
02003
02004 const TransitionProxy &tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
02005 if( ncs1[ipHi][ipLo] == 3 &&
02006 (fabs(tr.Emis().Aul()-1e-5) < 1e-8 ) )
02007 {
02008 ++nSkip;
02009 }
02010 else
02011 {
02012
02013
02014
02015 Save1LineData( tr , ioPUN , false , lgPrint);
02016 }
02017 }
02018 }
02019 fprintf( ioPUN , " %li lines skipped\n",nSkip);
02020 }
02021 }
02022
02023 return;
02024 }
02025
02026
02027
02028
02029 void FeIIPunDepart(
02030
02031 FILE* ioPUN ,
02032
02033 bool lgDoAll )
02034 {
02035
02036 # define NLEVDEP 11
02037
02038 const int LevDep[NLEVDEP]={1,10,25,45,64,124,206,249,295,347,371};
02039 long int i;
02040 static bool lgFIRST=true;
02041
02042 DEBUG_ENTRY( "FeIIPunDepart()" );
02043
02044
02045 if( lgFIRST && !lgDoAll )
02046 {
02047
02048 for( i=0; i<NLEVDEP; ++i )
02049 {
02050 fprintf( ioPUN , "%i\t", LevDep[i] );
02051 }
02052 fprintf( ioPUN , "\n");
02053 lgFIRST = false;
02054 }
02055
02056 if( lgDoAll )
02057 {
02058
02059 for( i=1; i<=FeII.nFeIILevel_local; ++i )
02060 {
02061 FeIIPun1Depart( ioPUN , i );
02062 fprintf( ioPUN , "\n");
02063 }
02064 }
02065 else
02066 {
02067
02068 for( i=0; i<NLEVDEP; ++i )
02069 {
02070 FeIIPun1Depart( ioPUN , LevDep[i] );
02071 fprintf( ioPUN , "\t");
02072 }
02073 fprintf( ioPUN , "\n");
02074 }
02075
02076 return;
02077 }
02078
02079
02080
02081
02082 void FeIIPun1Depart(
02083
02084 FILE * ioPUN ,
02085
02086 long int nPUN )
02087 {
02088
02089 DEBUG_ENTRY( "FeIIPun1Depart()" );
02090
02091 ASSERT( nPUN > 0 );
02092
02093 assert( ioPUN != NULL );
02094
02095
02096
02097 if( nPUN <= FeII.nFeIILevel_local )
02098 {
02099 fprintf( ioPUN , "%e ",Fe2DepCoef[nPUN-1] );
02100 }
02101 else
02102 {
02103 fprintf( ioPUN , "%e ",0. );
02104 }
02105
02106 return;
02107 }
02108
02109
02110
02111
02112 void FeIIReset(void)
02113 {
02114
02115 DEBUG_ENTRY( "FeIIReset()" );
02116
02117
02118 FeII.nFeIILevel_local = FeII.nFeIILevel_malloc;
02119
02120 return;
02121 }
02122
02123
02124
02125
02126
02127 void FeIIZero(void)
02128 {
02129
02130 DEBUG_ENTRY( "FeIIZero()" );
02131
02132
02133 FeII.Fe2_large_cool = 0.;
02134 FeII.ddT_Fe2_large_cool = 0.;
02135 FeII.Fe2_large_heat = 0.;
02136
02137
02138 FeII.lgLyaPumpOn = true;
02139
02140
02141
02142
02143
02144
02145 FeII.lgFeIILargeOn = false;
02146
02147
02148 FeII.fe2ener[0] = 0.;
02149 FeII.fe2ener[1] = 1e8;
02150
02151
02152
02153 FeII.ipRedisFcnResonance = ipCRD;
02154
02155 FeII.ipRedisFcnSubordinate = ipCRDW;
02156
02157
02158 FeII.fe2thresh = 0.;
02159
02160
02161
02162 FeII.lgSlow = false;
02163
02164
02165 FeII.lgPrint = false;
02166
02167
02168 FeII.lgSimulate = false;
02169
02170
02171
02172 if( lgFeIIMalloc )
02173 {
02174
02175 FeII.nFeIILevel_local = FeII.nFeIILevel_malloc;
02176 }
02177 else
02178 {
02179
02180
02181
02182 FeII.nFeIILevel_local = 16;
02183 }
02184
02185 return;
02186 }
02187
02188
02189 void FeIIPunPop(
02190
02191 FILE* ioPUN ,
02192
02193 bool lgPunchRange ,
02194
02195 long int ipRangeLo ,
02196
02197 long int ipRangeHi ,
02198
02199 bool lgPunchDensity )
02200 {
02201
02202 # define NLEVPOP 11
02203
02204 const int LevPop[NLEVPOP]={1,10,25,45,64,124,206,249,295,347,371};
02205 long int i;
02206 static bool lgFIRST=true;
02207 realnum denominator = 1.f;
02208
02209 DEBUG_ENTRY( "FeIIPunPop()" );
02210
02211
02212 if( !lgPunchDensity )
02213 denominator = SDIV( dense.xIonDense[ipIRON][1] );
02214
02215
02216
02217 if( lgFIRST && !lgPunchRange )
02218 {
02219
02220 for( i=0; i<NLEVPOP; ++i )
02221 {
02222
02223 fprintf( ioPUN , "%i\t", LevPop[i] );
02224 }
02225 fprintf( ioPUN , "\n");
02226 lgFIRST = false;
02227 }
02228
02229 if( lgPunchRange )
02230 {
02231
02232 ASSERT( ipRangeLo>=0 && ipRangeLo<ipRangeHi );
02233
02234
02235
02236 long nHigh = MIN2(ipRangeHi, FeII.nFeIILevel_local );
02237 for( i=ipRangeLo; i<nHigh; ++i )
02238 {
02239
02240 fprintf( ioPUN , "%.3e\t", Fe2LevelPop[i]/denominator );
02241 }
02242 fprintf( ioPUN , "\n");
02243 }
02244 else
02245 {
02246
02247
02248 for( i=0; i<NLEVPOP; ++i )
02249 {
02250 fprintf( ioPUN , "%.3e\t", Fe2LevelPop[LevPop[i]-1]/denominator );
02251 }
02252 fprintf( ioPUN , "\n");
02253 }
02254
02255 return;
02256 }
02257
02258
02259
02260
02261 #if 0
02262 void FeIIPun1Pop(
02263
02264 FILE * ioPUN ,
02265
02266 long int nPUN )
02267 {
02268 DEBUG_ENTRY( "FeIIPun1Pop()" );
02269
02270 ASSERT( nPUN > 0 );
02271
02272 assert( ioPUN != NULL );
02273
02274
02275
02276
02277 if( nPUN <= FeII.nFeIILevel_local )
02278 {
02279 fprintf( ioPUN , "%e ",Fe2LevelPop[nPUN-1] );
02280 }
02281 else
02282 {
02283 fprintf( ioPUN , "%e ",0. );
02284 }
02285
02286 return;
02287 }
02288 #endif
02289
02290
02291
02292
02293 STATIC int FeIIBandsCreate(
02294
02295
02296
02297 const char chFile[] )
02298 {
02299
02300 char chLine[FILENAME_PATH_LENGTH_2];
02301 const char* chFile1;
02302 FILE *ioDATA;
02303
02304 bool lgEOL;
02305 long int i,k;
02306
02307
02308
02309 static bool lgCalled=false;
02310
02311 DEBUG_ENTRY( "FeIIBandsCreate()" );
02312
02313
02314 if( lgCalled )
02315 {
02316
02317 return 0;
02318 }
02319 lgCalled = true;
02320
02321
02322 chFile1 = ( strlen(chFile) == 0 ) ? "FeII_bands.ini" : chFile;
02323
02324
02325 if( trace.lgTrace )
02326 {
02327 fprintf( ioQQQ, " FeIICreate opening %s:", chFile1 );
02328 }
02329
02330 ioDATA = open_data( chFile1, "r" );
02331
02332
02333 nFeIIBands = 0;
02334
02335
02336 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
02337 {
02338 fprintf( ioQQQ, " FeIICreate could not read first line of %s.\n", chFile1 );
02339 return 1;
02340 }
02341 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
02342 {
02343
02344
02345 if( chLine[0] != '#')
02346 ++nFeIIBands;
02347 }
02348
02349
02350 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
02351 {
02352 fprintf( ioQQQ, " FeIICreate could not rewind %s.\n", chFile1 );
02353 return 1;
02354 }
02355
02356 FeII_Bands = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(nFeIIBands) );
02357
02358
02359 for( i=0; i<nFeIIBands; ++i )
02360 {
02361 FeII_Bands[i] = (realnum *)MALLOC(sizeof(realnum)*(unsigned)(3) );
02362 }
02363
02364
02365 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
02366 {
02367 fprintf( ioQQQ, " FeIICreate could not read first line of %s.\n", chFile1 );
02368 return 1;
02369 }
02370 {
02371 i = 1;
02372 const long int iyr = 9, imo=6 , idy = 11;
02373 long iyrread, imoread , idyread;
02374 iyrread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
02375 imoread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
02376 idyread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
02377
02378 if(( iyrread != iyr ) ||
02379 ( imoread != imo ) ||
02380 ( idyread != idy ) )
02381 {
02382 fprintf( ioQQQ,
02383 " PROBLEM FeIIBandsCreate: the version of %s is not the current version.\n", chFile1 );
02384 fprintf( ioQQQ,
02385 " FeIIBandsCreate: I expected the magic numbers %li %li %li but found %li %li %li.\n",
02386 iyr, imo , idy ,iyrread, imoread , idyread );
02387 return 1;
02388 }
02389 }
02390
02391
02392 k = 0;
02393 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
02394 {
02395
02396
02397 if( chLine[0] != '#')
02398 {
02399 i = 1;
02400 FeII_Bands[k][0] = (realnum)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
02401 if( lgEOL )
02402 {
02403 fprintf( ioQQQ, " There should have been a number on this band line 1. Sorry.\n" );
02404 fprintf( ioQQQ, "string==%s==\n" ,chLine );
02405 return 1;
02406 }
02407 FeII_Bands[k][1] = (realnum)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
02408 if( lgEOL )
02409 {
02410 fprintf( ioQQQ, " There should have been a number on this band line 2. Sorry.\n" );
02411 fprintf( ioQQQ, "string==%s==\n" ,chLine );
02412 return 1;
02413 }
02414 FeII_Bands[k][2] = (realnum)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
02415 if( lgEOL )
02416 {
02417 fprintf( ioQQQ, " There should have been a number on this band line 3. Sorry.\n" );
02418 fprintf( ioQQQ, "string==%s==\n" ,chLine );
02419 return 1;
02420 }
02421
02422
02423 ++k;
02424 }
02425 }
02426
02427 for( i=0; i<nFeIIBands; ++i )
02428 {
02429
02430 if( FeII_Bands[i][0] <=0. || FeII_Bands[i][1] <=0. || FeII_Bands[i][2] <=0. )
02431 {
02432 fprintf( ioQQQ, " FeII band %li has none positive entry.\n",i );
02433 return 1;
02434 }
02435
02436 if( FeII_Bands[i][1] >= FeII_Bands[i][2] )
02437 {
02438 fprintf( ioQQQ, " FeII band %li has improper bounds.\n" ,i);
02439 return 1;
02440 }
02441
02442 }
02443
02444 fclose(ioDATA);
02445
02446
02447 return 0;
02448 }
02449
02450
02451
02452 void AssertFeIIDep( double *pred , double *BigError , double *StdDev )
02453 {
02454 long int n;
02455 double arg ,
02456 error ,
02457 sum2;
02458
02459 DEBUG_ENTRY( "AssertFeIIDep()" );
02460
02461 if( FeII.lgSimulate || !lgFeIIEverCalled )
02462 {
02463
02464 *pred = 0.;
02465 *BigError = 0.;
02466 *StdDev = 0.;
02467
02468 return;
02469 }
02470
02471
02472 ASSERT( FeII.nFeIILevel_local > 0 );
02473
02474
02475 *BigError = 0.;
02476 *pred = 0.;
02477 sum2 = 0;
02478 for( n=0; n<FeII.nFeIILevel_local; ++n )
02479 {
02480
02481 *pred += Fe2DepCoef[n];
02482
02483
02484 error = fabs( Fe2DepCoef[n] -1. );
02485
02486 *BigError = MAX2( *BigError , error );
02487
02488
02489 sum2 += POW2( Fe2DepCoef[n] );
02490 }
02491
02492
02493 arg = sum2 - POW2( *pred ) / (double)FeII.nFeIILevel_local;
02494 ASSERT( (arg >= 0.) );
02495 *StdDev = sqrt( arg / (double)(FeII.nFeIILevel_local - 1.) );
02496
02497
02498 *pred /= (double)(FeII.nFeIILevel_local);
02499
02500 return;
02501 }
02502
02503
02504 void FeII_OTS( void )
02505 {
02506 long int ipLo ,
02507 ipHi;
02508
02509 DEBUG_ENTRY( "FeII_OTS()" );
02510
02511 for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
02512 {
02513 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_local; ipHi++ )
02514 {
02515 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
02516
02517 if( tr.ipCont() < 1)
02518 continue;
02519
02520
02521 tr.Emis().ots() =
02522 tr.Emis().Aul()*
02523 (*tr.Hi()).Pop()*
02524 tr.Emis().Pdest();
02525
02526 ASSERT( tr.Emis().ots() >= 0. );
02527
02528
02529 RT_OTS_AddLine(tr.Emis().ots(),
02530 tr.ipCont() );
02531 }
02532 }
02533
02534 return;
02535 }
02536
02537
02538
02539
02540
02541
02542 void FeII_RT_Out(void)
02543 {
02544 long int ipLo , ipHi;
02545
02546 DEBUG_ENTRY( "FeIIRTOut()" );
02547
02548
02549 if( dense.xIonDense[ipIRON][1] > 0. )
02550 {
02551
02552 for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
02553 {
02554 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_local; ipHi++ )
02555 {
02556 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
02557
02558 if( tr.ipCont() < 1)
02559 continue;
02560
02561
02562 tr.outline_resonance();
02563
02564 }
02565 }
02566 }
02567
02568 return;
02569 }
02570
02571
02572
02573
02574 STATIC void FeIIContCreate(
02575
02576 double xLamLow ,
02577
02578 double xLamHigh ,
02579
02580 long int ncell )
02581 {
02582
02583 double step , wl1;
02584 long int i;
02585
02586
02587
02588 static bool lgCalled=false;
02589
02590 DEBUG_ENTRY( "FeIIContCreate()" );
02591
02592
02593 if( lgCalled )
02594 {
02595
02596 return;
02597 }
02598 lgCalled = true;
02599
02600
02601 nFeIIConBins = ncell;
02602
02603 FeII_Cont = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(nFeIIConBins+1) );
02604
02605
02606 for( i=0; i<nFeIIConBins+1; ++i )
02607 FeII_Cont[i] = (realnum *)MALLOC(sizeof(realnum)*(unsigned)(3) );
02608
02609 step = log10( xLamHigh/xLamLow)/ncell;
02610 wl1 = log10( xLamLow);
02611 for( i=0; i<nFeIIConBins+1; ++i)
02612
02613 FeII_Cont[i][0] = (realnum)pow(10. , ( wl1 + i*step ) );
02614
02615 return;
02616 }
02617
02618
02619 void ParseAtomFeII( Parser &p )
02620 {
02621 DEBUG_ENTRY( "ParseAtomFeII()" );
02622
02623
02624 FeII.lgFeIILargeOn = true;
02625
02626 if( lgFeIIMalloc )
02627 {
02628
02629 FeII.nFeIILevel_local = FeII.nFeIILevel_malloc;
02630 }
02631 else
02632 {
02633
02634 FeII.nFeIILevel_local = NFE2LEVN;
02635 }
02636
02637
02638
02639 if( p.nMatch("LEVE") )
02640 {
02641
02642 if( !lgFeIIMalloc )
02643 {
02644
02645 FeII.nFeIILevel_local = (long int)p.getNumberCheck("LEVEL");
02646 if( FeII.nFeIILevel_local <16 )
02647 {
02648 fprintf( ioQQQ, " This would be too few levels, must have at least 16.\n" );
02649 cdEXIT(EXIT_FAILURE);
02650 }
02651 else if( FeII.nFeIILevel_local > NFE2LEVN )
02652 {
02653 fprintf( ioQQQ, " This would be too many levels.\n" );
02654 cdEXIT(EXIT_FAILURE);
02655 }
02656 }
02657 }
02658
02659
02660 else if( p.nMatch("SLOW") )
02661 {
02662 FeII.lgSlow = true;
02663 }
02664
02665
02666 else if( p.nMatch("REDI") )
02667 {
02668 int ipRedis=0;
02669
02670
02671
02672
02673
02674 if( p.nMatch(" PRD") )
02675 {
02676 ipRedis = ipPRD;
02677 }
02678
02679 else if( p.nMatch(" CRD") )
02680 {
02681 ipRedis = ipCRD;
02682 }
02683
02684 else if( p.nMatch("CRDW") )
02685 {
02686 ipRedis = ipCRDW;
02687 }
02688
02689
02690 else if( !p.nMatch("SHOW") )
02691 {
02692 fprintf(ioQQQ," There should have been a second keyword on this command.\n");
02693 fprintf(ioQQQ," Options are _PRD, _CRD, CRDW (_ is space). Sorry.\n");
02694 cdEXIT(EXIT_FAILURE);
02695 }
02696
02697
02698 if( p.nMatch("RESO") )
02699 {
02700 FeII.ipRedisFcnResonance = ipRedis;
02701 }
02702
02703 else if( p.nMatch("SUBO") )
02704 {
02705 FeII.ipRedisFcnSubordinate = ipRedis;
02706 }
02707
02708 else if( p.nMatch("SHOW") )
02709 {
02710 fprintf(ioQQQ," FeII resonance lines are ");
02711 if( FeII.ipRedisFcnResonance ==ipCRDW )
02712 {
02713 fprintf(ioQQQ,"complete redistribution with wings\n");
02714 }
02715 else if( FeII.ipRedisFcnResonance ==ipCRD )
02716 {
02717 fprintf(ioQQQ,"complete redistribution with core only.\n");
02718 }
02719 else if( FeII.ipRedisFcnResonance ==ipPRD )
02720 {
02721 fprintf(ioQQQ,"partial redistribution.\n");
02722 }
02723 else
02724 {
02725 fprintf(ioQQQ," PROBLEM Impossible value for ipRedisFcnResonance.\n");
02726 TotalInsanity();
02727 }
02728
02729 fprintf(ioQQQ," FeII subordinate lines are ");
02730 if( FeII.ipRedisFcnSubordinate ==ipCRDW )
02731 {
02732 fprintf(ioQQQ,"complete redistribution with wings\n");
02733 }
02734 else if( FeII.ipRedisFcnSubordinate ==ipCRD )
02735 {
02736 fprintf(ioQQQ,"complete redistribution with core only.\n");
02737 }
02738 else if( FeII.ipRedisFcnSubordinate ==ipPRD )
02739 {
02740 fprintf(ioQQQ,"partial redistribution.\n");
02741 }
02742 else
02743 {
02744 fprintf(ioQQQ," PROBLEM Impossible value for ipRedisFcnSubordinate.\n");
02745 TotalInsanity();
02746 }
02747 }
02748 else
02749 {
02750 fprintf(ioQQQ," here should have been a second keyword on this command.\n");
02751 fprintf(ioQQQ," Options are RESONANCE, SUBORDINATE. Sorry.\n");
02752 cdEXIT(EXIT_FAILURE);
02753 }
02754 }
02755
02756
02757 else if( p.nMatch("TRAC") )
02758 {
02759 FeII.lgPrint = true;
02760 }
02761
02762
02763 else if( p.nMatch("SIMU") )
02764 {
02765
02766 FeII.lgSimulate = true;
02767 }
02768
02769
02770 else if( p.nMatch("CONT") )
02771 {
02772
02773
02774
02775 FeII.feconwlLo = (realnum)p.getNumberCheck("low wavelength");
02776 FeII.feconwlHi = (realnum)p.getNumberCheck("high wavelength");
02777 FeII.nfe2con = (long int)p.getNumberCheck("number of intervals");
02778
02779 if( FeII.feconwlLo<=0. || FeII.feconwlHi<=0. || FeII.nfe2con<= 0 )
02780 {
02781 fprintf(ioQQQ," there are three numbers on the FeII continuum command, start and end wavelengths, and number of intervals.\n");
02782 fprintf(ioQQQ," all three must be greater than zero, sorry.\n");
02783 cdEXIT(EXIT_FAILURE);
02784 }
02785
02786 if( FeII.feconwlLo>= FeII.feconwlHi )
02787 {
02788 fprintf(ioQQQ," there are three numbers on the FeII continuum command, start and end wavelengths, and number of intervals.\n");
02789 fprintf(ioQQQ," the second wavelength must be greater than the first, sorry.\n");
02790 cdEXIT(EXIT_FAILURE);
02791 }
02792 }
02793
02794
02795
02796 return;
02797 }
02798
02799 void PunFeII( FILE * io )
02800 {
02801 long int n, ipHi;
02802
02803 DEBUG_ENTRY( "PunFeII()" );
02804
02805 for( n=0; n<FeII.nFeIILevel_local-1; ++n)
02806 {
02807 for( ipHi=n+1; ipHi<FeII.nFeIILevel_local; ++ipHi )
02808 {
02809 const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][n]];
02810 if( tr.ipCont() > 0 )
02811 fprintf(io,"%li\t%li\t%.2e\n", n , ipHi ,
02812 tr.Emis().TauIn() );
02813 }
02814 }
02815
02816 return;
02817 }
02818
02819
02820 void FeIIPunchLineStuff( FILE * io , realnum xLimit , long index)
02821 {
02822 long int n, ipHi;
02823
02824 DEBUG_ENTRY( "FeIIPunchLineStuff()" );
02825
02826 for( n=0; n<FeII.nFeIILevel_local-1; ++n)
02827 {
02828 for( ipHi=n+1; ipHi<FeII.nFeIILevel_local; ++ipHi )
02829 {
02830 Save1Line( Fe2LevN[ipFe2LevN[ipHi][n]] , io , xLimit , index, GetDopplerWidth(dense.AtomicWeight[ipIRON]) );
02831 }
02832 }
02833
02834 return;
02835 }
02836
02837
02838 double FeIIRadPress(void)
02839 {
02840
02841
02842 realnum smallfloat=SMALLFLOAT*10.f;
02843 double press,
02844 RadPres1;
02845 # undef DEBUGFE
02846 # ifdef DEBUGFE
02847 double RadMax;
02848 long ipLoMax , ipHiMax;
02849 # endif
02850 long int ipLo, ipHi;
02851
02852 DEBUG_ENTRY( "FeIIRadPress()" );
02853
02854 press = 0.;
02855 if( !lgFeIIEverCalled )
02856 return 0.;
02857
02858 # ifdef DEBUGFE
02859 RadMax = 0;
02860 ipLoMax = -1;
02861 ipHiMax = -1;
02862 # endif
02863
02864 for( ipHi=1; ipHi<FeII.nFeIILevel_local; ++ipHi )
02865 {
02866 for( ipLo=0; ipLo<ipHi; ++ipLo)
02867 {
02868 const TransitionProxy &tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
02869
02870 if( tr.ipCont() > 0 &&
02871 (*tr.Hi()).Pop() > 1e-30 )
02872 {
02873 if( (*tr.Hi()).Pop() > smallfloat &&
02874 tr.Emis().PopOpc() > smallfloat )
02875 {
02876 RadPres1 = PressureRadiationLine( tr, GetDopplerWidth(dense.AtomicWeight[ipIRON]) );
02877
02878 # ifdef DEBUGFE
02879 if( RadPres1 > RadMax )
02880 {
02881 RadMax = RadPres1;
02882 ipLoMax = ipLo;
02883 ipHiMax = ipHi;
02884 }
02885 # endif
02886 press += RadPres1;
02887 }
02888 }
02889 }
02890 }
02891
02892 # ifdef DEBUGFE
02893
02894 if( iteration > 1 || nzone > 1558 )
02895 {
02896 fprintf(ioQQQ,"DEBUG FeIIRadPress finds P(FeII) = %.2e %.2e %li %li width %.2e\n",
02897 press ,
02898 RadMax ,
02899 ipLoMax ,
02900 ipHiMax ,
02901 RT_LineWidth(Fe2LevN[ipFe2LevN[9][0]],GetDopplerWidth(dense.AtomicWeight[ipIRON])) );
02902 DumpLine( Fe2LevN.begin()+ipFe2LevN[9][0] );
02903 }
02904 # endif
02905 # undef DEBUGFE
02906
02907 return press;
02908 }
02909
02910 #if 0
02911
02912 double FeII_InterEnergy(void)
02913 {
02914 double energy;
02915
02916 DEBUG_ENTRY( "FeII_InterEnergy()" );
02917
02918
02919
02920
02921 energy = 0.;
02922 for( long ipHi=1; ipHi<FeII.nFeIILevel_local; ++ipHi )
02923 {
02924 const TransitionList::iterator&tr = Fe2LevN.begin()+ipFe2LevN[ipHi][0];
02925 if( (*(*tr).Hi()).Pop() > 1e-30 )
02926 {
02927 energy +=
02928 (*(*tr).Hi()).Pop() * (*tr).EnergyErg;
02929 }
02930 }
02931
02932 return energy;
02933 }
02934 #endif
02935
02936
02937 #if defined(__HP_aCC)
02938 # pragma OPT_LEVEL 1
02939 #endif
02940
02941 STATIC void FeIILyaPump(void)
02942 {
02943
02944 long int ipLo ,
02945 ipHi;
02946 double EnerLyaProf2,
02947 EnerLyaProf3,
02948 EnergyWN,
02949 Gup_ov_Glo,
02950 PhotOccNum_at_nu,
02951 PumpRate,
02952 de,
02953 FeIILineWidthHz,
02954 taux;
02955
02956 DEBUG_ENTRY( "FeIILyaPump()" );
02957
02958
02959
02960
02961
02962
02963 if( FeII.lgLyaPumpOn )
02964 {
02965
02966
02967
02968
02969
02970
02971 de = 1.37194e-06*hydro.HLineWidth*2.0/3.0;
02972
02973 EnerLyaProf1 = 82259.0 - de*2.0;
02974 EnerLyaProf2 = 82259.0 - de;
02975 EnerLyaProf3 = 82259.0 + de;
02976 EnerLyaProf4 = 82259.0 + de*2.0;
02977
02978
02979 if( iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop() > SMALLFLOAT )
02980 {
02981
02982 PhotOccNumLyaCenter =
02983 MAX2(0.,1.0- iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pelec_esc() -
02984 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pesc())/
02985 (iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop()/iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop()*3. - 1.0);
02986 }
02987 else
02988 {
02989
02990 PhotOccNumLyaCenter = 0.;
02991 }
02992 }
02993 else
02994 {
02995 PhotOccNumLyaCenter = 0.;
02996 de = 0.;
02997 EnerLyaProf1 = 0.;
02998 EnerLyaProf2 = 0.;
02999 EnerLyaProf3 = 0.;
03000 EnerLyaProf4 = 0.;
03001 }
03002
03003
03004 if( FeII.lgLyaPumpOn )
03005 {
03006 for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
03007 {
03008 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_local; ipHi++ )
03009 {
03010 const TransitionProxy&tr=Fe2LevN[ipFe2LevN[ipHi][ipLo]];
03011
03012
03013 if( iteration == 1 )
03014 {
03015 taux = tr.Emis().TauIn();
03016 }
03017 else
03018 {
03019 taux = tr.Emis().TauTot();
03020 }
03021
03022
03023 Gup_ov_Glo = (*tr.Hi()).g()/(*tr.Lo()).g();
03024
03025
03026 EnergyWN = tr.EnergyWN();
03027
03028 if( EnergyWN >= EnerLyaProf1 && EnergyWN <= EnerLyaProf4 && taux > 1 )
03029 {
03030
03031
03032
03033
03034
03035
03036
03037
03038
03039
03040
03041
03042
03043
03044
03045 if( EnergyWN < EnerLyaProf2 )
03046 {
03047
03048 PhotOccNum_at_nu = PhotOccNumLyaCenter*(EnergyWN - EnerLyaProf1)/ de;
03049 }
03050 else if( EnergyWN < EnerLyaProf3 )
03051 {
03052
03053 PhotOccNum_at_nu = PhotOccNumLyaCenter;
03054 }
03055 else
03056 {
03057
03058 PhotOccNum_at_nu = PhotOccNumLyaCenter*(EnerLyaProf4 - EnergyWN)/de;
03059 }
03060
03061
03062
03063
03064
03066
03067 FeIILineWidthHz = 1.e8/(EnergyWN*299792.5)*sqrt(log(taux))*GetDopplerWidth(dense.AtomicWeight[ipIRON]);
03068
03069
03070 PumpRate = FeIILineWidthHz * PhotOccNum_at_nu * tr.Emis().Aul()*
03071 powi(82259.0f/EnergyWN,3);
03072
03073 PumpRate = tr.Emis().Aul()* PhotOccNum_at_nu;
03074
03075
03076 Fe2LPump[ipHi][ipLo] += PumpRate;
03077
03078
03079 Fe2LPump[ipLo][ipHi] += PumpRate*Gup_ov_Glo;
03080 }
03081 }
03082 }
03083 }
03084
03085 return;
03086 }
03087
03088
03089 #if defined(__HP_aCC)
03090 #pragma OPTIMIZE OFF
03091 #pragma OPTIMIZE ON
03092 #endif