26 #define INTBIG 2000000000
168 FILE *ioDATA, *ioROMAS , *ioBEHAR=NULL;
181 bool lgTooBig =
false;
189 fprintf(
ioQQQ,
" This is the second or later calculation in a grid.\n");
190 fprintf(
ioQQQ,
" The number of zones has been increased beyond what it was on the first calculation.\n");
191 fprintf(
ioQQQ,
" This can\'t be done since space has already been allocated.\n");
192 fprintf(
ioQQQ,
" Have the first calculation do the largest number of zones so that an increase is not needed.\n");
193 fprintf(
ioQQQ,
" Sorry.\n");
283 for( ipZ=0; ipZ<
LIMELM;++ipZ )
292 for( ipZ=0; ipZ< (LIMELM+3);++ipZ )
298 for( ion=0; ion < (LIMELM+1); ++ion )
342 fprintf(
ioQQQ,
" atmdat_readin opening level1.dat:");
349 fprintf(
ioQQQ,
" atmdat_readin could not read first line of level1.dat.\n");
356 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
360 if( chLine[0] !=
'#')
376 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
378 fprintf(
ioQQQ,
" atmdat_readin could not rewind level1.dat.\n");
385 fprintf(
ioQQQ,
" atmdat_readin could not read first line of level1.dat.\n");
396 if( ( nelem != 5 ) || ( nelec != 12 ) || ( ion != 19 ) )
399 " atmdat_readin: the version of level1.dat is not the current version.\n" );
401 " Please obtain the current version from the Cloudy web site.\n" );
403 " I expected to find the number 04 11 5 and got %li %li %li instead.\n" ,
404 nelem , nelec , ion );
405 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
412 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
414 if( i >= (nLevel1+1) )
416 fprintf(
ioQQQ,
" version conflict.\n");
419 if( chLine[0] !=
'#')
434 strncpy( chS2 , chLine , 2);
452 " atmdat_readin could not identify chem symbol on this level 1line:\n");
453 fprintf(
ioQQQ,
"%s\n", chLine );
454 fprintf(
ioQQQ,
"looking for this string==%2s==\n",chS2 );
466 fprintf(
ioQQQ,
" There should have been a number on this level1 line 1. Sorry.\n" );
467 fprintf(
ioQQQ,
"string==%s==\n" ,chLine );
477 fprintf(
ioQQQ,
" There should have been a number on this level1 line 2. Sorry.\n" );
478 fprintf(
ioQQQ,
"string==%s==\n" ,chLine );
485 fprintf(
ioQQQ,
" There should have been a number on this level1 line 3. Sorry.\n" );
486 fprintf(
ioQQQ,
"string==%s==\n" ,chLine );
493 fprintf(
ioQQQ,
" There should have been a number on this level1 line 4. Sorry.\n" );
494 fprintf(
ioQQQ,
"string==%s==\n" ,chLine );
501 fprintf(
ioQQQ,
" There should have been a number on this level1 line 5. Sorry.\n" );
502 fprintf(
ioQQQ,
"string==%s==\n" ,chLine );
511 fprintf(
ioQQQ,
" There should have been a number on this level1 line 6. Sorry.\n" );
512 fprintf(
ioQQQ,
"string==%s==\n" ,chLine );
523 fprintf(
ioQQQ,
" There should have been a number on this level1 line 7. Sorry.\n" );
524 fprintf(
ioQQQ,
"string==%s==\n" ,chLine );
534 fprintf(
ioQQQ,
" There should have been a number on this level1 line 8. Sorry.\n" );
535 fprintf(
ioQQQ,
"string==%s==\n" ,chLine );
583 fprintf(
ioQQQ,
" atmdat_readin opening level2.dat:");
590 fprintf(
ioQQQ,
" level2.dat error getting magic number\n" );
603 if( lgEOL || ( nelem != 6 ) || ( nelec != 8 ) || ( ion != 10 ) )
606 " atmdat_readin: the version of level2.dat is not the current version.\n" );
608 " I expected to find the number 06 08 10 and got %li %li %li instead.\n" ,
609 nelem , nelec , ion );
610 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
611 fprintf(
ioQQQ,
"Please obtain the correct version.\n");
617 while( i < nWindLine )
623 fprintf(
ioQQQ,
" level2.dat error getting line %li\n", i );
630 sscanf( chLine ,
"%lf %lf %lf %lf %lf %lf %lf " ,
653 fprintf(
ioQQQ,
" level2.dat error getting last magic number\n" );
656 sscanf( chLine ,
"%ld" , &magic2 );
659 fprintf(
ioQQQ,
" level2.dat ends will wrong magic number=%ld \n",
665 fprintf(
ioQQQ,
" OK \n");
689 fprintf(
ioQQQ,
" atmdat_readin opening UTA_Gu06.dat:");
691 ioGU06 =
open_data(
"UTA_Gu06.dat",
"r" );
694 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioGU06 ) != NULL )
698 if( chLine[0] !=
'#')
706 if( fseek( ioGU06 , 0 , SEEK_SET ) != 0 )
708 fprintf(
ioQQQ,
" atmdat_readin could not rewind UTA_Gu06.dat.\n");
715 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioGU06 ) != NULL )
719 if( chLine[0] !=
'#')
730 if( ( nelem != 2007 ) || ( nelec != 1 ) || ( ion != 23 ) )
733 " atmdat_readin: the version of UTA_Gu06.dat is not the current version.\n" );
735 " I expected to find the number 2007 1 23 and got %li %li %li instead.\n" ,
736 nelem , nelec , ion );
737 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
750 fprintf(
ioQQQ,
" atmdat_readin opening UTA_Behar.dat:");
752 ioBEHAR =
open_data(
"UTA_Behar.dat",
"r" );
756 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioBEHAR ) == NULL )
758 fprintf(
ioQQQ,
" atmdat_readin could not read first line of UTA_Behar.dat.\n");
772 if( ( nelem != 2002 ) || ( nelec != 10 ) || ( ion != 22 ) )
775 " atmdat_readin: the version of UTA_Behar.dat is not the current version.\n" );
777 " I expected to find the number 2002 10 22 and got %li %li %li instead.\n" ,
778 nelem , nelec , ion );
779 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
785 fprintf(
ioQQQ,
" atmdat_readin UTA data files open OK\n");
794 for( nelem=ipISO; nelem<
LIMELM; ++nelem )
804 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioBEHAR ) != NULL )
808 if( chLine[0] !=
'#')
812 if( fseek( ioBEHAR , 0 , SEEK_SET ) != 0 )
814 fprintf(
ioQQQ,
" atmdat_readin could not rewind UTA_Behar.dat.\n");
818 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioBEHAR ) == NULL )
820 fprintf(
ioQQQ,
" atmdat_readin could not read first line of UTA_Behar.dat.\n");
830 fprintf(
ioQQQ,
" atmdat_readin opening UTA_Kisielius.dat:");
832 ioROMAS =
open_data(
"UTA_Kisielius.dat",
"r" );
835 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioROMAS ) != NULL )
839 if( chLine[0] !=
'#')
849 for( i=0; i<nUTA_Behar+nUTA_Romas+nUTA_Gu06; i++ )
858 for( ion=0; ion<
ipIRON; ++ion )
860 nFeIonRomas[ion] = 0;
861 WLloRomas[ion] = 1e10;
864 nFeIonBehar[ion] = 0;
865 WLloBehar[ion] = 1e10;
876 if( fseek( ioROMAS , 0 , SEEK_SET ) != 0 )
878 fprintf(
ioQQQ,
" atmdat_readin could not rewind UTA_Kisielius.dat.\n");
887 for( nelem=ipISO; nelem<
LIMELM; ++nelem )
890 double sum_spon_auto;
893 double alam[
N]={ 1.859035 , 4.692138 , 10.324543 , 19.294364 ,
894 40.401292 , 44.442754 , 72.959719 };
895 double blam[
N]={-9.7855968,-11.159739, -12.914931 , -14.987272,
896 -18.853446 , -19.271781 ,-23.828388 };
897 double clam[
N]={ 8.2874628, 8.3043002, 8.3164038 , 8.3269937,
898 8.4312895 , 8.3730893 , 8.493802 };
901 double aA[
N] = {1.9067 , 1.8715 , 2.1033, 2.2815 ,
902 1.9511 , 1.9966 , 2.0852};
903 double bA[
N] = {8.4538 , 8.5528, 7.616, 6.9247 ,
904 7.9479, 7.9777 , 7.8349 };
906 double aDep[
N] = {29.54 , 12.07 , 24.23 , 7.35 , 7.37 , 11.14 , 11.99 };
907 double bDep[
N] = {-14.853 , -7.606 , -7.844 , 9.987 , 10.503 , 8.541 , 8.865 };
912 ion = nelem + 1 - ipISO;
920 f1 = alam[ipISO-2]*1e-4 + blam[ipISO-2]*1e-4*(nelem+1) +
921 clam[ipISO-2]*1e-4*
POW2(nelem+1.);
927 f1 = aA[ipISO-2]*log((
double)(nelem+1)) + bA[ipISO-2];
938 f1 = log((
double)(nelem+1));
940 sum_spon_auto = aDep[ipISO-2]*1e-2*f1 + bDep[ipISO-2]*1e-1;
944 sum_spon_auto = exp( sum_spon_auto );
945 sum_spon_auto = pow(10., sum_spon_auto )*1e13;
949 sum_spon_auto = pow(10., sum_spon_auto )*1e13;
963 fprintf(
ioQQQ,
" %2s%2li\t%.4f\t%.3e\t%.3e\t%.3e\n",
982 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
985 if( chLine[0] !=
'#')
988 double WLAng, Aul, oscill , Aauto;
990 sscanf( chLine ,
"%4li%5li%8lf%13le%12le",
991 &ion,&i2,&WLAng,&Aul,&Aauto);
992 sscanf( &chLine[54] ,
"%13le", &oscill);
998 ASSERT( ion >= 0 && ion < ipIRON );
1014 WLloGu[ion-1] =
MIN2( WLAng , WLloGu[ion-1] );
1015 WLhiGu[ion-1] =
MAX2( WLAng , WLloGu[ion-1] );
1029 double frac_ioniz = Aauto/(Aul + Aauto);
1030 ASSERT( frac_ioniz >=0. && frac_ioniz <=1. );
1047 enum {DEBUG_LOC=
false};
1051 fprintf(
ioQQQ,
"DEBUG Gu UTA\t%li\t%.3f\t%.3e\t%.3e\t%.3e\n",
1072 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
1075 if( chLine[0] !=
'#')
1077 long int i1, i2, i3;
1078 double f1, f2, oscill;
1081 sscanf( chLine ,
"%li\t%li\t%li\t%lf\t%le\t%le\t%le",
1082 &i1,&i2,&i3,&f1,&f2,&frac_relax,&oscill );
1088 ASSERT( i1 >= 0 && i1 < ipIRON );
1112 WLloBehar[i1] =
MIN2( f1 , WLloBehar[i1] );
1113 WLhiBehar[i1] =
MAX2( f1 , WLhiBehar[i1] );
1132 ASSERT( frac_relax >=0. && frac_relax <=1. );
1135 fprintf(
ioQQQ,
"data set %li line %li %2s%2li\t%.4f\t%.3e\t%.3e\n",
1157 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
1160 if( chLine[0] !=
'#')
1162 long int i1, i2, i3;
1163 double f1, f2, oscill;
1166 sscanf( chLine ,
"%li\t%li\t%li\t%lf\t%le\t%le\t%le",
1167 &i1,&i2,&i3,&f1,&f2,&frac_relax,&oscill );
1173 ASSERT( i1 >= 0 && i1 < ipIRON );
1197 WLloRomas[i1] =
MIN2( f1 , WLloRomas[i1] );
1198 WLhiRomas[i1] =
MAX2( f1 , WLhiRomas[i1] );
1216 ASSERT( frac_relax >=0. && frac_relax <=1. );
1225 ASSERT( i == (nUTA_Behar+nUTA_Romas+nUTA_Gu06) );
1229 fprintf(
ioQQQ ,
"summary of UTA data sets read in\nion numb WLlo WLhi\n");
1230 fprintf(
ioQQQ ,
"Behar 01 total lines %li\n" , nUTA_Behar );
1231 for( ion=0; ion<
ipIRON;++ion )
1233 if( nFeIonBehar[ion] )
1235 fprintf(
ioQQQ,
"%3li %6li %7.3f %7.3f \n",
1243 fprintf(
ioQQQ ,
"Gu 06 total lines %li\n" , nUTA_Gu06 );
1244 for( ion=0; ion<
ipIRON;++ion )
1248 fprintf(
ioQQQ,
"%3li %6li %7.3f %7.3f \n",
1256 fprintf(
ioQQQ ,
"Romas Kisielius 03 total lines %li\n", nUTA_Romas );
1257 for( ion=0; ion<
ipIRON;++ion )
1259 if( nFeIonRomas[ion] )
1261 fprintf(
ioQQQ,
"%3li %6li %7.3f %7.3f \n",
1278 for( i=0; i<nUTA_Behar+nUTA_Gu06; ++i )
1280 if(
UTALines[i].Hi->nelem-1 == ipIRON )
1287 if( nFeIonRomas[
UTALines[i].Hi->IonStg-1] )
1294 nUTA = nUTA_Behar + nUTA_Gu06 + nUTA_Romas;
1301 nUTA = nUTA_Behar + nUTA_Gu06;
1352 fprintf(
ioQQQ,
" atmdat_readin opening mewe_gbar.dat:");
1354 ioDATA =
open_data(
"mewe_gbar.dat",
"r" );
1357 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) == NULL )
1359 fprintf(
ioQQQ,
" mewe_gbar.dat error getting magic number\n" );
1363 sscanf( chLine ,
"%ld" , &magic1 );
1364 if( magic1 != 9101 )
1366 fprintf(
ioQQQ,
" mewe_gbar.dat starts with wrong magic number=%ld \n",
1372 for( i=1; i < 210; i++ )
1374 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) == NULL )
1376 fprintf(
ioQQQ,
" mewe_gbar.dat error getting line %li\n", i );
1381 sscanf( chLine,
"%lf %lf %lf %lf ", &help[0], &help[1], &help[2], &help[3] );
1382 for(
int l=0; l < 4; ++l )
1387 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) == NULL )
1389 fprintf(
ioQQQ,
" mewe_gbar.dat error getting last magic number\n" );
1393 sscanf( chLine ,
"%ld" , &magic2 );
1395 if( magic1 != magic2 )
1397 fprintf(
ioQQQ,
" mewe_gbar.dat ends will wrong magic number=%ld \n",
1405 fprintf(
ioQQQ,
" OK \n");
1410 for( nelem=0; nelem <
LIMELM; nelem++ )
1411 for( ion=0; ion <
LIMELM; ion++ )
1418 for( nelem=2; nelem < LIMELM; nelem++ )
1421 for( ion=0; ion <= nelem; ion++ )
1424 nelec = nelem - ion + 1;
1434 ASSERT( imax > 0 && imax <= 10 );
1460 for( ipZ=0; ipZ<
HS_NZ; ++ipZ )
1463 if( ipZ>1 && ipZ<5 )
continue;
1465 for( iCase=0; iCase<2; ++iCase )
1472 sprintf( chFilename,
"HS_e%ld%c.dat", ipZ+1, ( iCase == 0 ) ?
'a' :
'b' );
1475 fprintf(
ioQQQ,
" atmdat_readin opening Hummer Storey emission file %s:",chFilename);
1482 fprintf(
ioQQQ,
" OK\n");
1486 i = fscanf( ioDATA,
"%li %li ",
1495 for( ipTemp=0; ipTemp <
atmdat.
ntemp[iCase][ipZ]; ipTemp++ )
1499 long int junk, junk2 , ne;
1500 fscanf( ioDATA,
" %lf %li %lf %c %li %ld ",
1507 for( j=0; j < ne; j++ )
1509 fscanf( ioDATA,
"%lf ", &
atmdat.
Emiss[iCase][ipZ][ipTemp][ipDens][j] );
1522 fprintf(
ioQQQ,
"\n");
1523 for( ipTemp=0; ipTemp<
atmdat.
ntemp[iCase][ipZ]; ipTemp++ )
1527 fprintf(
ioQQQ,
"\n");
1550 for(
int ns=0; ns < 7; ns++ )
1553 for(
int nelec=0; nelec < 10; nelec++ )
1567 const char* chFilename;
1596 chFilename =
"mewe_nelectron.dat";
1599 fprintf(
ioQQQ,
" init_yield opening %s:", chFilename );
1619 while( ch1 ==
'#' || ch1 ==
'*' )
1623 fprintf(
ioQQQ,
" %s error getting line %i\n", chFilename, ns );
1629 sscanf( chLine,
"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
1630 &temp[0], &temp[1], &temp[2], &temp[3], &temp[4],
1631 &temp[5], &temp[6], &temp[7], &temp[8], &temp[9],
1632 &temp[10],&temp[11],&temp[12],&temp[13],&temp[14] );
1638 for(
int j=0; j < 10; j++ )
1655 fprintf(
ioQQQ,
" Auger yields will be killed.\n");
1656 fprintf(
ioQQQ,
" OK\n");
1661 chFilename =
"mewe_fluor.dat";
1664 fprintf(
ioQQQ,
" init_yield opening %s:", chFilename );
1676 fprintf(
ioQQQ,
" %s error getting line %i\n", chFilename, 0 );
1680 while( chLine[0] ==
'#' );
1688 int nDima[NKM] = { 0, 1, 2, 2, 3, 4, 4, 5, 5, 6 };
1695 sscanf( chLine,
"%lf %lf %lf %lf %lf %lf %lf",
1696 &temp[0], &temp[1], &temp[2], &temp[3], &temp[4],
1697 &temp[5], &temp[6] );
1714 nAuger = (int)temp[3];
1739 while( chLine[0]==
'#' && !lgEOL );