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