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