00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #include "cddefines.h"
00020 #include "cddrive.h"
00021 #include "physconst.h"
00022 #include "mean.h"
00023 #include "taulines.h"
00024 #include "struc.h"
00025 #include "iso.h"
00026 #include "mole.h"
00027 #include "hyperfine.h"
00028 #include "rt.h"
00029 #include "lines_service.h"
00030 #include "doppvel.h"
00031 #include "dense.h"
00032 #include "h2.h"
00033 #include "magnetic.h"
00034 #include "hydrogenic.h"
00035 #include "secondaries.h"
00036 #include "grainvar.h"
00037 #include "lines.h"
00038 #include "dynamics.h"
00039 #include "colden.h"
00040 #include "continuum.h"
00041 #include "ionbal.h"
00042 #include "yield.h"
00043 #include "prt.h"
00044 #include "iterations.h"
00045 #include "heavy.h"
00046 #include "conv.h"
00047 #include "geometry.h"
00048 #include "called.h"
00049 #include "helike.h"
00050 #include "opacity.h"
00051 #include "rfield.h"
00052 #include "phycon.h"
00053 #include "timesc.h"
00054 #include "radius.h"
00055 #include "atomfeii.h"
00056 #include "monitor_results.h"
00057 #include "thermal.h"
00058 #include "wind.h"
00059 #include "hmi.h"
00060 #include "pressure.h"
00061 #include "elementnames.h"
00062 #include "ipoint.h"
00063 #include "gammas.h"
00064 #include "atmdat.h"
00065 #include "hcmap.h"
00066 #include "input.h"
00067 #include "save.h"
00068 #include "optimize.h"
00069 #include "warnings.h"
00070 #include "grid.h"
00071 #include "mole_priv.h"
00072 #include "atomfeii.h"
00073
00074
00075 extern "C" int wavelength_compare (const void * a, const void * b)
00076 {
00077 LinSv a1 = *(LinSv*)a;
00078 LinSv b1 = *(LinSv*)b;
00079
00080 if( b1.wavelength > a1.wavelength )
00081 return 1;
00082 else if( b1.wavelength < a1.wavelength )
00083 return -1;
00084 else
00085 return 0;
00086 }
00087
00088
00089 STATIC void FindStrongestLineLabels( void )
00090 {
00091 long low_index=0;
00092 long high_index=0;
00093 long j_min = 0;
00094 double MaxFlux = 0.;
00095 long ipMaxFlux = 0;
00096 long j = 0;
00097
00098 ASSERT( LineSave.ipass==1 );
00099
00100 memcpy(LineSvSortWL, LineSv, (unsigned)LineSave.nsum*sizeof( LinSv ));
00101 qsort((void *)LineSvSortWL,(size_t)LineSave.nsum,sizeof( LinSv ),wavelength_compare);
00102
00103 while( rfield.anu[j_min]+0.5*rfield.widflx[j_min] < RYDLAM/LineSvSortWL[0].wavelength)
00104 j_min++;
00105
00106 for( j=0; j<rfield.nflux; j++ )
00107 {
00108 if( j < j_min )
00109 {
00110 strcpy( rfield.chLineLabel[j], " " );
00111 continue;
00112 }
00113
00114 ASSERT( LineSvSortWL[low_index].wavelength != 0. );
00115
00116 while( RYDLAM/LineSvSortWL[low_index].wavelength < rfield.anu[j]-0.5*rfield.widflx[j] && low_index < LineSave.nsum-1 )
00117 {
00118 low_index++;
00119 if( LineSvSortWL[low_index].wavelength == 0. )
00120 {
00121
00122 for( long j1=j; j1<rfield.nflux; j1++ )
00123 strcpy( rfield.chLineLabel[j1], " " );
00124 return;
00125 }
00126 }
00127
00128 high_index = low_index;
00129 ASSERT( LineSvSortWL[high_index].wavelength != 0. );
00130
00131 while( RYDLAM/LineSvSortWL[high_index].wavelength < rfield.anu[j]+0.5*rfield.widflx[j] && high_index < LineSave.nsum-1 )
00132 {
00133 high_index++;
00134 if( LineSvSortWL[high_index].wavelength == 0. )
00135 {
00136 high_index--;
00137 break;
00138 }
00139 }
00140
00141 high_index--;
00142
00143 ASSERT( LineSvSortWL[low_index].wavelength > 0. );
00144 ASSERT( LineSvSortWL[high_index].wavelength > 0. );
00145 ASSERT( RYDLAM/LineSvSortWL[low_index].wavelength >= rfield.anu[j]-0.5*rfield.widflx[j] );
00146 ASSERT( RYDLAM/LineSvSortWL[high_index].wavelength <= rfield.anu[j]+0.5*rfield.widflx[j] );
00147
00148 MaxFlux = 0.;
00149 ipMaxFlux = 0;
00150
00151 for( long k = low_index; k <= high_index; k++ )
00152 {
00153 if( strcmp( LineSvSortWL[k].chALab, "Coll" )==0 ||
00154 strcmp( LineSvSortWL[k].chALab, "Heat" )==0 ||
00155 strcmp( LineSvSortWL[k].chALab, "Pump" )==0 ||
00156 strcmp( LineSvSortWL[k].chALab, "pump" )==0 ||
00157 strcmp( LineSvSortWL[k].chALab, "nInu" )==0 ||
00158 strcmp( LineSvSortWL[k].chALab, "nFnu" )==0 ||
00159 strcmp( LineSvSortWL[k].chALab, "InwT" )==0 ||
00160 strcmp( LineSvSortWL[k].chALab, "InwC" )==0 ||
00161 strcmp( LineSvSortWL[k].chALab, "Inwd" )==0 ||
00162 strcmp( LineSvSortWL[k].chALab, "Ca A" )==0 ||
00163 strcmp( LineSvSortWL[k].chALab, "Ca B" )==0 ||
00164 strcmp( LineSvSortWL[k].chALab, "Pho+" )==0 ||
00165 strcmp( LineSvSortWL[k].chALab, "Pcon" )==0 ||
00166 strcmp( LineSvSortWL[k].chALab, "Q(H)" )==0 ||
00167 strcmp( LineSvSortWL[k].chALab, "Unit" )==0 )
00168 continue;
00169
00170 if( LineSvSortWL[k].SumLine[0] > MaxFlux )
00171 {
00172 MaxFlux = LineSvSortWL[k].SumLine[0];
00173 ipMaxFlux = k;
00174 }
00175 }
00176
00177
00178 if( ipMaxFlux > 0 )
00179 strcpy( rfield.chLineLabel[j], LineSvSortWL[ipMaxFlux].chALab );
00180 }
00181
00182 return;
00183 }
00184
00185 # ifdef USE_NLTE7
00186
00187 static void runNLTE(int ipPun)
00188 {
00189 long nelem = ipNEON;
00190 int nouter[11] = {0,2,3,3,3,4,3,5,5,6,0};
00191 ASSERT(dense.gas_phase[nelem] != 0.);
00192
00193 fprintf( save.ipPnunit[ipPun], "data Ferland University of Kentucky Cloudy 10 11-5-2011\n");
00194 fprintf( save.ipPnunit[ipPun], "case NeXY6\n");
00195 fprintf( save.ipPnunit[ipPun], "code Cloudy\n");
00196 fprintf( save.ipPnunit[ipPun], "atom Ne 10\n");
00197 fprintf( save.ipPnunit[ipPun], "calctime %11.4e %11.4e\n",4.5,30.);
00198
00200
00201 double avgchrg = 0.;
00202 double m2 = 0.;
00203 double m3 = 0.;
00204
00205 for( int ion = 1; ion < nelem-1; ion++)
00206 {
00207 avgchrg += (dense.xIonDense[nelem][ion]/dense.gas_phase[nelem] * ion);
00208 }
00209
00210 for( int ion = 1; ion < nelem-1; ion++)
00211 {
00212 m2 += dense.xIonDense[nelem][ion]/dense.gas_phase[nelem]*POW2((ion - avgchrg));
00213 m3 += dense.xIonDense[nelem][ion]/dense.gas_phase[nelem]*POW3((ion - avgchrg));
00214
00215 }
00216
00217 double Eint = 0.;
00218 double partfun = 0.;
00219 long count_nrglvls = 0;
00220 for( int ion = 1; ion < nelem-1; ion++)
00221 {
00222
00223 long ipSpec = 0;
00224 for( int i = 0; i < nSpecies; i++)
00225 {
00226 ipSpec = i;
00227 if(dBaseStates[i][0].nelem()-1 == nelem && dBaseStates[i][0].IonStg()-1 == ion)
00228 break;
00229 ipSpec = -1;
00230 }
00231
00232 if( ipSpec != -1)
00233 {
00234 for( int lvl = 0; lvl < dBaseSpecies[ipSpec].numLevels_max; lvl++)
00235 {
00236
00237 Eint += dBaseStates[ipSpec][lvl].Pop()*dBaseStates[ipSpec][lvl].energy().eV();
00238 partfun += dBaseStates[ipSpec][lvl].g()*exp(-1*dBaseStates[ipSpec][lvl].energy().eV()/
00239 phycon.te_eV);
00240 count_nrglvls += 1;
00241 }
00242 }
00243 else
00244 {
00245 printf("\nNot using :\t%li\t%i\n",nelem+1,ion+1);
00246
00247 }
00248 }
00249 Eint = Eint / dense.gas_phase[nelem];
00250
00251
00252 fprintf( save.ipPnunit[ipPun],"\nsummary_quantities\n");
00253 fprintf( save.ipPnunit[ipPun],"plasma %11.4e %11.4e\n",phycon.te_eV,dense.eden);
00254 fprintf( save.ipPnunit[ipPun],"time %11.4e\n",0.);
00255 fprintf( save.ipPnunit[ipPun],"zbar %11.4e\n",avgchrg);
00256 fprintf( save.ipPnunit[ipPun],"m2 %11.4e\n",m2);
00257 fprintf( save.ipPnunit[ipPun],"m3 %11.4e\n",m3);
00258
00259 fprintf( save.ipPnunit[ipPun],"eint %11.4e\n",0.);
00260 fprintf( save.ipPnunit[ipPun],"deintdt %11.4e\n",0.);
00261 fprintf( save.ipPnunit[ipPun],"pfn %11.4e\n",partfun);
00262 fprintf( save.ipPnunit[ipPun],"nmax_eff %i\n",0);
00263 fprintf( save.ipPnunit[ipPun],"ploss %11.4e %11.4e %11.4e %9.3e\n",0.,0.,0.,thermal.totcol);
00264
00265 fprintf( save.ipPnunit[ipPun],"\nion_stages %li\n",nelem-2);
00266
00267
00268 double PIR[nelem+1];
00269
00270 double autoion[nelem+1];
00271
00272 for( int ion = 1; ion < nelem-1; ion++)
00273 {
00274 double tempRecomb = 0.;
00275 double f_aauto = 0.;
00276 double f_aphoto = 0.;
00277 double f_acoll = 0.;
00278 double f_Scoll = 0.;
00279 double f_Sphoto = 0.;
00280 double f_auto = 0.;
00281 autoion[ion] = 0.;
00282 PIR[ion] = 0.;
00283
00284 if( ion != 0 && ionbal.RateRecomTot[nelem][ion-1] != 0.)
00285 {
00286 tempRecomb = ionbal.RateRecomTot[nelem][ion-1];
00287 f_aauto = dense.eden*ionbal.DR_Badnell_rate_coef[nelem][ion-1]/tempRecomb;
00288 f_aphoto = dense.eden*ionbal.RR_rate_coef_used[nelem][ion-1]/tempRecomb;
00289 f_acoll = dense.eden*ionbal.CotaRate[ion-1]/tempRecomb;
00290 }
00291
00292 for( int i=0; i < Heavy.nsShells[nelem][ion]; i++)
00293 {
00294
00295 if( ion != nelem+1)
00296 {
00297 PIR[ion] += ionbal.PhotoRate_Shell[nelem][ion][i][0];
00298 }
00299 }
00300
00301 if( ion != nelem+1 && ionbal.RateIonizTot(nelem,ion) != 0.)
00302 {
00303 f_Scoll = ionbal.CollIonRate_Ground[nelem][ion][0]/ionbal.RateIonizTot(nelem,ion);
00304 f_Sphoto = PIR[ion]/ionbal.RateIonizTot(nelem,ion);
00305 f_auto = autoion[ion]/ionbal.RateIonizTot(nelem,ion);
00306 }
00307
00308 autoion[ion] = ionbal.UTA_ionize_rate[nelem][ion] + secondaries.csupra[nelem][ion];
00309 if( ion > 0)
00310 {
00311 autoion[ion] += 0.0;
00312 }
00313 fprintf( save.ipPnunit[ipPun],"ion %li %11.4e %i\n"
00314 " %11.3e %11.3e %11.3e %11.3e\n"
00315 " %11.3e %11.3e %11.3e %11.3e\n\n",
00316 nelem+1-ion,
00317 dense.xIonDense[nelem][ion]/dense.gas_phase[nelem],
00318 nouter[ion],
00319 ionbal.RateIonizTot(nelem,ion),
00320 f_Scoll,
00321 f_Sphoto,
00322 f_auto,
00323 tempRecomb,
00324 f_acoll,
00325 f_aphoto,
00326 f_aauto);
00327 }
00328
00329
00330 fprintf( save.ipPnunit[ipPun],"\nenergy_levels %li\n",count_nrglvls);
00331
00332 double GcollBF = 0.0;
00333
00334 double GphotoBF = 0.0;
00335
00336 double GautoBF = 0.0;
00337
00338 double Gtotal = 0.0;
00339
00340 double QcollBF = 0.0;
00341
00342 double QphotoBF = 0.0;
00343
00344 double QautoBF = 0.0;
00345
00346 double Qtotal = 0.0;
00347
00348 double totallevelpop = 0.;
00349
00350 for( int ion = 1; ion < nelem-1; ion++)
00351 {
00352
00353 long ipSpec = 0;
00354 for( int i = 0; i < nSpecies; i++)
00355 {
00356 ipSpec = i;
00357 if(dBaseStates[i][0].nelem()-1 == nelem && dBaseStates[i][0].IonStg()-1 == ion)
00358 break;
00359 ipSpec = -1;
00360 }
00361 if( ipSpec != -1)
00362 {
00363 for( int lvl = 0; lvl < dBaseSpecies[ipSpec].numLevels_max; lvl++)
00364 {
00365
00366 totallevelpop += dBaseStates[ipSpec][lvl].Pop();
00367 }
00368 }
00369 }
00370
00371 for( int ion = 1; ion < nelem-1; ion++)
00372 {
00373
00374 long ipSpec = 0;
00375 for( int i = 0; i < nSpecies; i++)
00376 {
00377 ipSpec = i;
00378 if(dBaseStates[i][0].nelem()-1 == nelem && dBaseStates[i][0].IonStg()-1 == ion)
00379 break;
00380 ipSpec = -1;
00381 }
00382
00383
00384 double ELCF = 0.;
00385 for( int i = 0; i < ion; i++ )
00386 {
00387 ELCF += Heavy.Valence_IP_Ryd[nelem][i]*EVRYD;
00388 }
00389
00390 if( ipSpec != -1)
00391 {
00392
00393
00394 for( int lvl = 0; lvl < dBaseSpecies[ipSpec].numLevels_max; lvl++)
00395 {
00396
00397
00398 if( lvl == 0 && ion != nelem+1)
00399 {
00400 GcollBF = ionbal.CollIonRate_Ground[nelem][ion][0];
00401 GphotoBF = PIR[ion];
00402 GautoBF = autoion[ion];
00403 }
00404 else
00405 {
00406 GcollBF = 0.;
00407 GphotoBF = 0.;
00408 GautoBF = 0.;
00409 }
00410
00411 Gtotal = dBaseStates[ipSpec][lvl].DestCollBB() +
00412 dBaseStates[ipSpec][lvl].DestPhotoBB() +
00413 GcollBF + GphotoBF + GautoBF;
00414
00415 double f_GcollBB,
00416 f_GphotoBB,
00417 f_GcollBF,
00418 f_GphotoBF,
00419 f_GautoBF;
00420
00421 if( Gtotal != 0.)
00422 {
00423 f_GcollBB = dBaseStates[ipSpec][lvl].DestCollBB()/Gtotal;
00424 f_GphotoBB = dBaseStates[ipSpec][lvl].DestPhotoBB()/Gtotal;
00425 f_GcollBF = GcollBF/Gtotal;
00426 f_GphotoBF = GphotoBF/Gtotal;
00427 f_GautoBF = GautoBF/Gtotal;
00428 }
00429 else
00430 {
00431 f_GcollBB = 0.;
00432 f_GphotoBB = 0.;
00433 f_GcollBF = 0.;
00434 f_GphotoBF = 0.;
00435 f_GautoBF = 0.;
00436 }
00437
00438 if( lvl == 0 && ion != 0 )
00439 {
00440 QautoBF = dense.eden*ionbal.DR_Badnell_rate_coef[nelem][ion-1];
00441 QphotoBF = dense.eden*ionbal.RR_rate_coef_used[nelem][ion-1];
00442 QcollBF = dense.eden*ionbal.CotaRate[ion-1];
00443 }
00444 else
00445 {
00446 QautoBF = 0.;
00447 QphotoBF = 0.;
00448 QcollBF = 0.;
00449 }
00450
00451
00452 Qtotal = dBaseStates[ipSpec][lvl].CreatCollBB() +
00453 0.0 + QcollBF + QphotoBF + QautoBF;
00454
00455 double f_QcollBB,
00456 f_QphotoBB,
00457 f_QcollBF,
00458 f_QphotoBF,
00459 f_QautoBF;
00460
00461 if( Qtotal != 0.)
00462 {
00463 f_QcollBB = dBaseStates[ipSpec][lvl].CreatCollBB()/Qtotal;
00464 f_QphotoBB = 0.0/Qtotal;
00465 f_QcollBF = QcollBF/Qtotal;
00466 f_QphotoBF = QphotoBF/Qtotal;
00467 f_QautoBF = QautoBF/Qtotal;
00468 }
00469 else
00470 {
00471 f_QcollBB = 0.;
00472 f_QphotoBB = 0.;
00473 f_QcollBF = 0.;
00474 f_QphotoBF = 0.;
00475 f_QautoBF = 0.;
00476 }
00477
00478
00479 fprintf( save.ipPnunit[ipPun],"elev %li %3i %11.4e %11.6e %11.4e\n"
00480 " %11.3e %11.3e %11.3e %11.3e %11.3e %11.3e\n"
00481 " %11.3e %11.3e %11.3e %11.3e %11.3e %11.3e\n"
00482 " %i %i %i\n\n",
00483 nelem+1-ion,
00484 lvl+1,
00485 dBaseStates[ipSpec][lvl].g(),
00486 dBaseStates[ipSpec][lvl].energy().eV()+ELCF,
00487 dBaseStates[ipSpec][lvl].Pop()/totallevelpop,
00488 Gtotal*dBaseStates[ipSpec][lvl].Pop(),
00489 f_GcollBB,
00490 f_GphotoBB,
00491 f_GcollBF,
00492 f_GphotoBF,
00493 f_GautoBF,
00494 Qtotal,
00495 f_QcollBB,
00496 f_QphotoBB,
00497 f_QcollBF,
00498 f_QphotoBF,
00499 f_QautoBF,
00500 0,
00501 0,
00502 nouter[ion]);
00503 }
00504 }
00505 else
00506 printf("\nNot using :\t%li\t%i for levels\n",nelem+1,ion+1);
00507 }
00508 return;
00509 }
00510 # endif
00511
00512
00513
00514
00515 inline realnum PrettyTransmission(long j, realnum transmission)
00516 {
00517 if( save.ResolutionAbs < realnum(0.) )
00518
00519 return transmission;
00520 else
00521 {
00522 realnum corr = save.ResolutionAbs*rfield.widflx[j]/rfield.anu[j];
00523 return realnum(max(0., 1. - (1.-transmission)*corr ));
00524 }
00525 }
00526
00527
00528
00529 STATIC void SaveResults1Line(
00530 FILE * ioPUN,
00531
00532 const char *chLab,
00533 realnum wl,
00534 double xInten,
00535 const char *chFunction);
00536
00537
00538 STATIC void SaveGaunts(FILE* ioPUN);
00539
00540
00541
00542 STATIC void SaveResults(FILE* ioPUN);
00543
00544 STATIC void SaveLineStuff(
00545 FILE * ioPUN,
00546 const char *chJob ,
00547 realnum xLimit);
00548
00549
00550 STATIC void AGN_Hemis(FILE *ioPUN );
00551
00552
00553 STATIC void SaveNewContinuum(FILE * ioPUN );
00554
00555
00556 STATIC void SaveLineIntensity(FILE * ioPUN , long int ipPun, realnum Threshold);
00557
00558 char *chDummy;
00559
00560
00561 STATIC realnum SaveFeII_cont( long int ipCont , long ipFeII_Cont_type )
00562 {
00563
00564 DEBUG_ENTRY( "SaveFeII_cont()" );
00565
00566
00567 if( ipFeII_Cont_type == 3 )
00568 return( FeII_Cont[ipCont][1]+FeII_Cont[ipCont][2] );
00569 else
00570 return(FeII_Cont[ipCont][ipFeII_Cont_type]);
00571 }
00572
00573 void SaveDo(
00574
00575 const char *chTime)
00576 {
00577 bool lgTlkSav;
00578 long int
00579 ipPun,
00580 i,
00581 i1,
00582 ion,
00583 ipConMax,
00584 ipHi,
00585 ipLinMax,
00586 ipLo,
00587 ips,
00588 j,
00589 jj,
00590 limit,
00591 nelem;
00592 double ConMax,
00593 RateInter,
00594 av,
00595 conem,
00596 eps,
00597 flxatt,
00598 flxcor,
00599 flxin,
00600 flxref,
00601 flxtrn,
00602 fout,
00603 fref,
00604 fsum,
00605 opConSum,
00606 opLinSum,
00607 stage,
00608 sum,
00609 texc,
00610 xLinMax;
00611
00612 DEBUG_ENTRY( "SaveDo()" );
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636 if( save.nsave < 1 )
00637 {
00638 return;
00639 }
00640
00641
00642
00643 if( grid.lgGrid )
00644 {
00645 if( chTime[0]=='L' )
00646 GridGatherInCloudy();
00647 }
00648
00649
00650
00651
00652
00653 if( strcmp(chTime,"LAST") == 0 )
00654 {
00655
00656 FindStrongestLineLabels();
00657 }
00658
00659 for( ipPun=0; ipPun < save.nsave; ipPun++ )
00660 {
00661 if( save.lgPunHeader[ipPun] && !nMatch(save.chHeader[ipPun],save.chNONSENSE) )
00662 {
00663 fprintf( save.ipPnunit[ipPun], "%s", save.chHeader[ipPun] );
00664 save.lgPunHeader[ipPun] = false;
00665 }
00666
00667
00668 save.ipConPun = ipPun;
00669
00670
00671 bool lgNoHitFirstBranch = false;
00672
00673
00674
00675
00676
00677
00678
00679 if( iterations.lgLastIt || !save.lgPunLstIter[ipPun] ||
00680
00681
00682
00683
00684 ( lgAbort && strcmp(chTime,"LAST") == 0 ) ||
00685 ( dynamics.lgTimeDependentStatic && dynamics.lgStatic_completed ) )
00686 {
00687
00688 if( strcmp(save.chSave[ipPun],"ABUN") == 0 )
00689 {
00690
00691 if( strcmp(chTime,"LAST") != 0 )
00692 {
00693 fprintf( save.ipPnunit[ipPun], "%.2f",
00694 log10(MAX2(SMALLFLOAT,dense.gas_phase[ipHYDROGEN])) );
00695 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
00696 {
00697
00698
00699 fprintf( save.ipPnunit[ipPun], "\t%.2f",
00700 log10(MAX2(SMALLFLOAT,dense.gas_phase[nelem])) );
00701 }
00702 fprintf( save.ipPnunit[ipPun], "\n" );
00703 }
00704 }
00705
00706 else if( strcmp(save.chSave[ipPun],"21CM") == 0 )
00707 {
00708
00709 if( strcmp(chTime,"LAST") != 0 )
00710 {
00711 fprintf( save.ipPnunit[ipPun],
00712 "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
00713
00714 radius.depth_mid_zone,
00715 hyperfine.Tspin21cm ,
00716 phycon.te ,
00717
00718 3.84e-7* iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauCon() /
00719 SDIV( HFLines[0].Emis().TauCon() ),
00720
00721 (*HFLines[0].Lo()).Pop() ,
00722 (*HFLines[0].Hi()).Pop() ,
00723 OccupationNumberLine( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) ),
00724 HFLines[0].Emis().TauCon() ,
00725 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauCon(),
00726 HFLines[0].Emis().PopOpc(),
00727
00728 (iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop())/SDIV( hyperfine.Tspin21cm),
00729
00730
00731 HFLines[0].Emis().TauIn(),
00732 TexcLine( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) ) ,
00733 colden.H0_ov_Tspin,
00734
00735 colden.H0_21cm_lower,
00736 colden.H0_21cm_upper,
00737 -0.068/log((colden.H0_21cm_upper/3.)/colden.H0_21cm_lower)
00738 );
00739 }
00740 }
00741
00742 else if( strcmp(save.chSave[ipPun],"AGES") == 0 )
00743 {
00744
00745 if( strcmp(chTime,"LAST") != 0 )
00746 {
00747 int ipCO, ipOH;
00748 ipCO = findspecies("CO")->index;
00749 ipOH = findspecies("OH")->index;
00750 fprintf( save.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00751
00752 radius.depth_mid_zone,
00753
00754 dense.pden*BOLTZMANN*1.5*phycon.te/ thermal.htot,
00755
00756 timesc.time_H2_Dest_here,
00757
00758 1./SDIV((ipCO != -1) ? mole.species[ipCO].snk : 0.),
00759
00760 1./SDIV((ipOH != -1) ? mole.species[ipOH].snk : 0.),
00761
00762 1./(dense.eden*2.90e-10/(phycon.te70*phycon.te10/phycon.te03)) );
00763 }
00764 }
00765
00766 else if( strcmp(save.chSave[ipPun]," AGN") == 0 )
00767 {
00768 if( strcmp(chTime,"LAST") == 0 )
00769 {
00770 if( strcmp( save.chSaveArgs[ipPun], "HECS" ) == 0 )
00771 {
00772
00773 AGN_He1_CS(save.ipPnunit[ipPun]);
00774 }
00775 if( strcmp( save.chSaveArgs[ipPun], "HEMI" ) == 0 )
00776 {
00777
00778 AGN_Hemis(save.ipPnunit[ipPun]);
00779 }
00780 else
00781 {
00782 fprintf( ioQQQ, " SaveDo does not recognize flag %4.4s set for AGN save. This is impossible.\n",
00783 save.chSave[ipPun] );
00784 ShowMe();
00785 cdEXIT(EXIT_FAILURE);
00786 }
00787 }
00788 }
00789
00790 else if( strcmp(save.chSave[ipPun],"MONI") == 0 )
00791 {
00792 if( strcmp(chTime,"LAST") == 0 )
00793 {
00794
00795 lgCheckMonitors(save.ipPnunit[ipPun]);
00796 }
00797 }
00798
00799 else if( strcmp(save.chSave[ipPun],"AVER") == 0 )
00800 {
00801 if( strcmp(chTime,"LAST") == 0 )
00802 {
00803
00804 save_average( ipPun );
00805 }
00806 }
00807
00808 else if( strncmp(save.chSave[ipPun],"CHA",3) == 0 )
00809 {
00810 if( strcmp(chTime,"LAST") == 0 )
00811 {
00812
00813 ChargTranPun( save.ipPnunit[ipPun] , save.chSave[ipPun] );
00814 }
00815 }
00816
00817 else if( strcmp( save.chSave[ipPun],"CHIA") == 0)
00818 {
00819
00820 int ipLo = 0;
00821 int ipHi = 0;
00822 double fupsilon = 0.;
00823 double initTemp = 4.0;
00824 double finalTemp = 9.1;
00825 double stepTemp = 0.2;
00826 fprintf( save.ipPnunit[ipPun],"Species\tLo\tHi\tWLAng");
00827 for(double logtemp = initTemp;logtemp < finalTemp;logtemp = logtemp + stepTemp )
00828 {
00829 fprintf( save.ipPnunit[ipPun],"\t%2.1f",logtemp);
00830 }
00831 fprintf( save.ipPnunit[ipPun],"\n");
00832 for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
00833 {
00834 for( EmissionList::iterator tr=dBaseTrans[ipSpecies].Emis().begin();
00835 tr != dBaseTrans[ipSpecies].Emis().end(); ++tr)
00836 {
00837 if( dBaseTrans[ipSpecies].chLabel() == "Chianti" )
00838 {
00839 ipLo = tr->Tran().ipLo();
00840 ipHi = tr->Tran().ipHi();
00841 fprintf( save.ipPnunit[ipPun],"%s\t%i\t%i\t",
00842 dBaseSpecies[ipSpecies].chLabel,ipLo+1,ipHi+1);
00843 fprintf( save.ipPnunit[ipPun],"%5.3e",tr->Tran().WLAng());
00844 for(double logtemp = initTemp;logtemp < finalTemp;logtemp = logtemp + stepTemp )
00845 {
00846 fupsilon = CHIANTI_Upsilon(ipSpecies, ipELECTRON, ipHi, ipLo, pow(10.,logtemp));
00847 fprintf( save.ipPnunit[ipPun],"\t%5.3e",fupsilon);
00848 }
00849 fprintf( save.ipPnunit[ipPun],"\n");
00850 }
00851 }
00852 }
00853 }
00854
00855 else if( strcmp(save.chSave[ipPun],"COOL") == 0 ||
00856 strcmp(save.chSave[ipPun],"EACH") == 0 )
00857 {
00858
00859 if( strcmp(chTime,"LAST") != 0 )
00860 CoolSave(save.ipPnunit[ipPun], save.chSave[ipPun]);
00861 }
00862
00863 else if( strcmp(save.chSave[ipPun],"DOMI") == 0 )
00864 {
00865
00866 if( strcmp(chTime,"LAST") != 0 )
00867 {
00868 molecule *debug_species = findspecies( save.chSpeciesDominantRates[ipPun] );
00869 if (debug_species == null_mole)
00870 {
00871 fprintf( ioQQQ,"Error in SAVE DOMINANT RATES, species %s not found\n",
00872 save.chSpeciesDominantRates[ipPun]);
00873 }
00874 else
00875 {
00876 fprintf( save.ipPnunit[ipPun],
00877 "%e\t%e\t", radius.depth_mid_zone, mole.species[ debug_species->index ].column );
00878 mole_dominant_rates( debug_species, save.ipPnunit[ipPun] );
00879 }
00880 }
00881 }
00882
00883 else if( strncmp(save.chSave[ipPun],"DYN" , 3) == 0 )
00884 {
00885
00886 if( strcmp(chTime,"LAST") != 0 )
00887 DynaSave(save.ipPnunit[ipPun] ,save.chSave[ipPun][3] );
00888 }
00889
00890 else if( strcmp(save.chSave[ipPun],"ENTH") == 0 )
00891 {
00892 if( strcmp(chTime,"LAST") != 0 )
00893 fprintf( save.ipPnunit[ipPun],
00894 "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00895 radius.depth_mid_zone,
00896 phycon.EnthalpyDensity,
00897 phycon.EnergyExcitation,
00898 phycon.EnergyIonization,
00899 phycon.EnergyBinding ,
00900 0.5*POW2(wind.windv)*dense.xMassDensity ,
00901 5./2.*pressure.PresGasCurr ,
00902 magnetic.EnthalpyDensity);
00903 }
00904
00905 else if( strcmp(save.chSave[ipPun],"COLU") == 0 )
00906 {
00907
00908 if( strcmp(chTime,"LAST") == 0 )
00909 {
00910 PrtColumns(save.ipPnunit[ipPun],"TABLE",ipPun);
00911 }
00912 }
00913
00914 else if( strcmp(save.chSave[ipPun],"COLS") == 0 )
00915 {
00916
00917 if( strcmp(chTime,"LAST") == 0 )
00918 {
00919
00920 save_colden( save.ipPnunit[ipPun] );
00921 }
00922 }
00923
00924 else if( strcmp(save.chSave[ipPun],"COMP") == 0 )
00925 {
00926
00927 if( strcmp(chTime,"LAST") != 0 )
00928 {
00929 for( jj=0; jj<rfield.nflux; jj = jj + save.ncSaveSkip)
00930 {
00931 fprintf( save.ipPnunit[ipPun], "%10.2e%10.2e%10.2e\n",
00932 rfield.anu[jj], rfield.comup[jj]/rfield.widflx[jj],
00933 rfield.comdn[jj]/rfield.widflx[jj] );
00934 }
00935 }
00936 }
00937
00938
00939 else if( strcmp(save.chSave[ipPun],"CON ") == 0 )
00940 {
00941
00942
00943
00944
00945 bool lgPrintThis =false;
00946 if( save.lgSaveEveryZone[ipPun] )
00947 {
00948
00949 if( strcmp(chTime,"LAST") != 0 )
00950 {
00951
00952 if( nzone==1 )
00953 {
00954 lgPrintThis = true;
00955 }
00956 else if( nzone%save.nSaveEveryZone[ipPun]==0 )
00957 {
00958 lgPrintThis = true;
00959 }
00960 }
00961 else
00962 {
00963
00964 if( nzone!=1 && nzone%save.nSaveEveryZone[ipPun]!=0 )
00965 {
00966 lgPrintThis = true;
00967 }
00968 }
00969 }
00970 else
00971 {
00972
00973 if( strcmp(chTime,"LAST") == 0 )
00974 lgPrintThis = true;
00975 }
00976 ASSERT( !save.lgSaveEveryZone[ipPun] || save.nSaveEveryZone[ipPun]>0 );
00977 if( lgPrintThis )
00978 {
00979 if( save.lgSaveEveryZone[ipPun] && nzone!=1)
00980 fprintf( save.ipPnunit[ipPun], "%s\n",
00981 save.chHashString );
00982
00983
00984 int nEmType = (int)save.punarg[ipPun][0];
00985 ASSERT( nEmType==0 || nEmType==1 );
00986
00987 const realnum *trans_coef_total=rfield.getCoarseTransCoef();
00988 for( j=0; j<rfield.nflux; j = j+save.ncSaveSkip)
00989 {
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000 flxin = ((double)rfield.flux_total_incident[nEmType][j])
01001 *rfield.anu2[j]*
01002 EN1RYD/rfield.widflx[j];
01003
01004
01005 realnum resolution = ( save.Resolution < realnum(0.) ) ?
01006 rfield.anu[j]/rfield.widflx[j] : save.Resolution;
01007
01008
01009 flxref = (rfield.anu2[j]*((double)rfield.ConRefIncid[nEmType][j]+rfield.ConEmitReflec[nEmType][j])/rfield.widflx[j] +
01010 rfield.anu[j]*resolution*rfield.reflin[nEmType][j])*EN1RYD;
01011
01012
01013 flxatt = ((double)rfield.flux[nEmType][j])*rfield.anu2[j]*EN1RYD/
01014 rfield.widflx[j]*radius.r1r0sq *
01015 PrettyTransmission( j, trans_coef_total[j] );
01016
01017
01018 conem = ((double)rfield.ConEmitOut[nEmType][j]/
01019 rfield.widflx[j]*rfield.anu2[j] + resolution*
01020 rfield.outlin[nEmType][j]*rfield.anu[j])*radius.r1r0sq*
01021 EN1RYD*geometry.covgeo;
01022
01023
01024 flxtrn = conem + flxatt;
01025
01026
01027 fprintf( save.ipPnunit[ipPun],"%.5e\t", AnuUnit(rfield.AnuOrg[j]) );
01028
01029 fprintf( save.ipPnunit[ipPun],"%.3e\t", flxin );
01030
01031 fprintf( save.ipPnunit[ipPun],"%.3e\t", flxatt );
01032
01033 fprintf( save.ipPnunit[ipPun],"%.3e\t", conem );
01034
01035 fprintf( save.ipPnunit[ipPun],"%.3e\t", flxtrn );
01036
01037 fprintf( save.ipPnunit[ipPun],"%.3e\t", flxref );
01038
01039 fprintf( save.ipPnunit[ipPun],"%.3e\t", flxref + flxtrn );
01040
01041
01042 fprintf( save.ipPnunit[ipPun],"%.3e\t", rfield.anu[j]*resolution*rfield.reflin[nEmType][j]*EN1RYD );
01043
01044
01045 fprintf( save.ipPnunit[ipPun],"%.3e\t", rfield.anu[j]*resolution*rfield.outlin[nEmType][j]*EN1RYD*
01046 radius.r1r0sq*geometry.covgeo );
01047
01048 fprintf( save.ipPnunit[ipPun], "%4.4s\t%4.4s\t",
01049
01050 rfield.chLineLabel[j] ,
01051
01052 rfield.chContLabel[j] );
01053
01054
01055 fprintf( save.ipPnunit[ipPun], "%.2f\n", rfield.line_count[j]/rfield.widflx[j]*rfield.anu[j] );
01056 }
01057 }
01058 }
01059
01060
01061
01062 else if( strcmp(save.chSave[ipPun],"CONN") == 0 )
01063 {
01064 if( strcmp(chTime,"LAST") == 0 )
01065 SaveNewContinuum( save.ipPnunit[ipPun]);
01066 }
01067
01068 else if( strcmp(save.chSave[ipPun],"CONC") == 0 )
01069 {
01070
01071
01072
01073 if( strcmp(chTime,"LAST") == 0 )
01074 {
01075
01076 for( j=0; j<rfield.nflux; j = j + save.ncSaveSkip)
01077 {
01078 flxin = rfield.flux_total_incident[0][j]*rfield.anu2[j]*
01079 EN1RYD/rfield.widflx[j];
01080
01081 fprintf( save.ipPnunit[ipPun], "%.5e\t%.3e\n",
01082 AnuUnit(rfield.AnuOrg[j]), flxin );
01083 }
01084 }
01085 }
01086
01087 else if( strcmp(save.chSave[ipPun],"CONG") == 0 )
01088 {
01089
01090 if( strcmp(chTime,"LAST") == 0 )
01091 {
01092
01093
01094 for( j=0; j<rfield.nflux; j = j + save.ncSaveSkip)
01095 {
01096 fsum = gv.GraphiteEmission[j]*rfield.anu2[j]*
01097 EN1RYD/rfield.widflx[j] *radius.r1r0sq*geometry.covgeo;
01098
01099 fout = gv.SilicateEmission[j]*rfield.anu2[j]*
01100 EN1RYD/rfield.widflx[j] *radius.r1r0sq*geometry.covgeo;
01101
01102
01103
01104 fprintf( save.ipPnunit[ipPun], "%.5e\t%.3e\t%.3e\t%.3e\n",
01105 AnuUnit(rfield.AnuOrg[j]) , fsum , fout ,
01106
01107
01108 gv.GrainEmission[j]*rfield.anu2[j]*
01109 EN1RYD/rfield.widflx[j] *radius.r1r0sq*geometry.covgeo );
01110 }
01111 }
01112 }
01113
01114 else if( strcmp(save.chSave[ipPun],"CONR") == 0 )
01115 {
01116
01117
01118
01119 if( strcmp(chTime,"LAST") == 0 )
01120 {
01121 if( geometry.lgSphere )
01122 {
01123 fprintf( save.ipPnunit[ipPun], " Reflected continuum not predicted when SPHERE is set.\n" );
01124 fprintf( ioQQQ ,
01125 "\n\n>>>>>>>>>>>>>\n Reflected continuum not predicted when SPHERE is set.\n" );
01126 cdEXIT(EXIT_FAILURE);
01127 }
01128
01129 for( j=0; j<rfield.nflux; j = j + save.ncSaveSkip)
01130 {
01131
01132 realnum resolution = ( save.Resolution < realnum(0.) ) ?
01133 rfield.anu[j]/rfield.widflx[j] : save.Resolution;
01134
01135
01136 flxref = rfield.anu2[j]*((double)rfield.ConRefIncid[0][j]+rfield.ConEmitReflec[0][j])/
01137 rfield.widflx[j]*EN1RYD;
01138
01139 fref = rfield.anu[j]*resolution*rfield.reflin[0][j]*EN1RYD;
01140
01141 if( rfield.flux_total_incident[0][j] > 1e-25 )
01142 {
01143 av = rfield.ConRefIncid[0][j]/rfield.flux_total_incident[0][j];
01144 }
01145 else
01146 {
01147 av = 0.;
01148 }
01149
01150 fprintf( save.ipPnunit[ipPun], "%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4s\n",
01151 AnuUnit(rfield.AnuOrg[j]), flxref, fref, flxref + fref,
01152 av, rfield.chContLabel[j] );
01153 }
01154 }
01155 }
01156
01157 else if( strcmp(save.chSave[ipPun],"CNVE") == 0 )
01158 {
01159
01160 if( strcmp(chTime,"LAST") != 0 )
01161 {
01162 fprintf( save.ipPnunit[ipPun],
01163 "%.5e\t%li\t%.4e\t%.4f\t%.4e\t%.4e\t%.3f\t%.4e\t%.4e\t%.4f\n",
01164 radius.depth_mid_zone,
01165 conv.nPres2Ioniz,
01166 pressure.PresTotlCurr,
01167 pressure.PresTotlError*100.,
01168 dense.EdenTrue,
01169 dense.eden,
01170 (dense.EdenTrue - dense.eden)*100./dense.EdenTrue,
01171 thermal.htot,
01172 thermal.ctot,
01173 (thermal.htot - thermal.ctot)*100./thermal.htot );
01174 }
01175 }
01176
01177 else if( strcmp(save.chSave[ipPun],"CONB") == 0 )
01178 {
01179
01180
01181
01182 if( strcmp(chTime,"LAST") != 0 )
01183 {
01184 for( j=0; j<rfield.nupper; j = j + save.ncSaveSkip)
01185 {
01186 fprintf( save.ipPnunit[ipPun], "%14.5e %14.5e %14.5e\n",
01187 AnuUnit(rfield.AnuOrg[j]), rfield.anu[j], rfield.widflx[j] );
01188 }
01189 }
01190 }
01191
01192 else if( strcmp(save.chSave[ipPun],"COND") == 0 )
01193 {
01194 ASSERT( fp_equal( save.punarg[ipPun][0] , (realnum)2. ) ||
01195 fp_equal( save.punarg[ipPun][0] , (realnum)1.) );
01196
01197
01198 if( strcmp(chTime,"LAST") == 0 &&
01199 fp_equal(save.punarg[ipPun][0] , (realnum)1.) )
01200 {
01201
01202
01203 for( j=0; j<rfield.nflux; j = j+save.ncSaveSkip)
01204 {
01205
01206 realnum resolution = ( save.Resolution < realnum(0.) ) ?
01207 rfield.anu[j]/rfield.widflx[j] : save.Resolution;
01208
01209 double EmisLin = resolution*EN1RYD*
01210 rfield.DiffuseLineEmission[j]*rfield.anu[j];
01211 double EmisCon = rfield.ConEmitLocal[nzone][j]*
01212 rfield.anu2[j]*EN1RYD/rfield.widflx[j];
01213 fprintf( save.ipPnunit[ipPun], "%.5e\t%.5e\t%.5e\t%.5e\n",
01214 AnuUnit(rfield.AnuOrg[j]),
01215 EmisCon ,
01216 EmisLin ,
01217 EmisCon+EmisLin);
01218 }
01219 }
01220 else if( strcmp(chTime,"LAST") != 0 &&
01221 fp_equal(save.punarg[ipPun][0] , (realnum)2.) )
01222 {
01223
01224 static bool lgMustPrintHeader=true;
01225 if( lgMustPrintHeader )
01226 {
01227 lgMustPrintHeader = false;
01228 fprintf( save.ipPnunit[ipPun], "%.5e",
01229 AnuUnit(rfield.AnuOrg[0]) );
01230 for( j=1; j<rfield.nflux; j = j+save.ncSaveSkip)
01231 {
01232 fprintf( save.ipPnunit[ipPun], "\t%.5e",
01233 AnuUnit(rfield.AnuOrg[j]) );
01234 }
01235 fprintf( save.ipPnunit[ipPun], "\n" );
01236 }
01237
01238 realnum resolution = ( save.Resolution < realnum(0.) ) ?
01239 rfield.anu[0]/rfield.widflx[0] : save.Resolution;
01240 double EmisLin = resolution*EN1RYD*
01241 rfield.DiffuseLineEmission[0]*rfield.anu[0];
01242 double EmisCon = rfield.ConEmitLocal[nzone][0]*
01243 rfield.anu2[0]*EN1RYD/rfield.widflx[0];
01244 fprintf( save.ipPnunit[ipPun], "%.5e",
01245 EmisCon+EmisLin);
01246 for( j=1; j<rfield.nflux; j = j+save.ncSaveSkip)
01247 {
01248
01249 resolution = ( save.Resolution < realnum(0.) ) ?
01250 rfield.anu[j]/rfield.widflx[j] : save.Resolution;
01251 double EmisLin = resolution*EN1RYD*
01252 rfield.DiffuseLineEmission[j]*rfield.anu[j];
01253 double EmisCon = rfield.ConEmitLocal[nzone][j]*
01254 rfield.anu2[j]*EN1RYD/rfield.widflx[j];
01255 fprintf( save.ipPnunit[ipPun], "\t%.5e",
01256 EmisCon+EmisLin);
01257 }
01258 fprintf( save.ipPnunit[ipPun], "\n" );
01259 }
01260 }
01261
01262 else if( strcmp(save.chSave[ipPun],"CONE") == 0 )
01263 {
01264
01265
01266
01267 if( strcmp(chTime,"LAST") == 0 )
01268 {
01269
01270 for( j=0; j<rfield.nflux;j = j +save.ncSaveSkip)
01271 {
01272
01273 realnum resolution = ( save.Resolution < realnum(0.) ) ?
01274 rfield.anu[j]/rfield.widflx[j] : save.Resolution;
01275
01276
01277 flxref = (rfield.anu2[j]*(rfield.ConRefIncid[0][j]+rfield.ConEmitReflec[0][j])/
01278 rfield.widflx[j] + rfield.anu[j]*resolution*
01279 rfield.reflin[0][j])*EN1RYD;
01280
01281
01282 conem = (rfield.ConEmitOut[0][j])/rfield.widflx[j]*rfield.anu2[j] +
01283 resolution*rfield.outlin[0][j]*rfield.anu[j];
01284
01285 conem *= radius.r1r0sq*EN1RYD*geometry.covgeo;
01286
01287
01288
01289 fprintf( save.ipPnunit[ipPun], "%.5e\t%.3e\t%.3e\t%.3e\t%4.4s\t%4.4s\n",
01290 AnuUnit(rfield.AnuOrg[j]),
01291 flxref,
01292 conem,
01293 flxref + conem,
01294 rfield.chLineLabel[j],
01295 rfield.chContLabel[j]
01296 );
01297 }
01298 }
01299 }
01300
01301
01302 else if( strcmp(save.chSave[ipPun],"CONf") == 0 )
01303 {
01304 if( strcmp(chTime,"LAST") == 0 )
01305 {
01306 long nu_hi , nskip;
01307 if( save.punarg[ipPun][0] > 0. )
01308
01309
01310 j = ipFineCont( save.punarg[ipPun][0] );
01311 else
01312 j = 0;
01313
01314
01315 if( save.punarg[ipPun][1]> 0. )
01316 nu_hi = ipFineCont( save.punarg[ipPun][1]);
01317 else
01318 nu_hi = rfield.nfine;
01319
01320
01321 nskip = (long)save.punarg[ipPun][2];
01322
01323 do
01324 {
01325 realnum sum1 = rfield.fine_opt_depth[j];
01326 realnum xnu = rfield.fine_anu[j];
01327 for( jj=1; jj<nskip; ++jj )
01328 {
01329 xnu += rfield.fine_anu[j+jj];
01330 sum1 += rfield.fine_opt_depth[j+jj];
01331 }
01332 fprintf( save.ipPnunit[ipPun],
01333 "%.6e\t%.3e\n",
01334 AnuUnit(xnu/nskip),
01335 sexp(sum1/nskip) );
01336 j += nskip;
01337 } while( j < nu_hi );
01338 }
01339 }
01340
01341 else if( strcmp(save.chSave[ipPun],"CONi") == 0 )
01342 {
01343
01344
01345
01346
01347
01348 if( strcmp(chTime,"LAST") != 0 )
01349 {
01350
01351 if( save.punarg[ipPun][0] <= 0. )
01352 {
01353 i1 = 1;
01354 }
01355 else if( save.punarg[ipPun][0] < 100. )
01356 {
01357 i1 = ipoint(save.punarg[ipPun][0]);
01358 }
01359 else
01360 {
01361 i1 = (long int)save.punarg[ipPun][0];
01362 }
01363
01364 fref = 0.;
01365 fout = 0.;
01366 fsum = 0.;
01367 sum = 0.;
01368 flxin = 0.;
01369
01370 for( j=i1-1; j < rfield.nflux; j++ )
01371 {
01372 fref += rfield.flux[0][j]*opac.opacity_abs[j];
01373 fout += rfield.otslin[j]*opac.opacity_abs[j];
01374 fsum += rfield.otscon[j]*opac.opacity_abs[j];
01375 sum += rfield.ConInterOut[j]*opac.opacity_abs[j];
01376 flxin += (rfield.outlin[0][j] + rfield.outlin_noplot[j])*opac.opacity_abs[j];
01377 }
01378 fprintf( save.ipPnunit[ipPun], "%10.2e%10.2e%10.2e%10.2e%10.2e\n",
01379 fref, fout, fsum, sum, flxin );
01380 }
01381 }
01382
01383 else if( strcmp(save.chSave[ipPun],"CONI") == 0 )
01384 {
01385
01386
01387
01388
01389 if( save.lgSaveEveryZone[ipPun] || (strcmp(chTime,"LAST") == 0) )
01390 {
01391
01392 bool lgPrt=false;
01393 if( save.lgSaveEveryZone[ipPun] )
01394 fprintf(save.ipPnunit[ipPun],"#save every zone %li\n", nzone);
01395
01396
01397
01398
01399 if( save.punarg[ipPun][0] <= 0. )
01400 i1 = 1;
01401 else if( save.punarg[ipPun][0] < 100. )
01402 i1 = ipoint(save.punarg[ipPun][0]);
01403 else
01404 i1 = (long int)save.punarg[ipPun][0];
01405
01406 sum = 0.;
01407 for( j=i1-1; j < rfield.nflux; j++ )
01408 {
01409 flxcor = rfield.flux[0][j] +
01410 rfield.otslin[j] +
01411 rfield.otscon[j] +
01412 rfield.ConInterOut[j] +
01413 rfield.outlin[0][j] + rfield.outlin_noplot[j];
01414
01415 sum += flxcor*opac.opacity_abs[j];
01416 }
01417
01418 if( sum > 0. )
01419 sum = 1./sum;
01420 else
01421 sum = 1.;
01422
01423 fsum = 0.;
01424
01425 for( j=i1-1; j<rfield.nflux; ++j)
01426 {
01427 flxcor = rfield.flux[0][j] +
01428 rfield.otslin[j] +
01429 rfield.otscon[j] +
01430 rfield.ConInterOut[j]+
01431 rfield.outlin[0][j] + rfield.outlin_noplot[j];
01432
01433 fsum += flxcor*opac.opacity_abs[j];
01434
01435
01436
01437 RateInter = flxcor*opac.opacity_abs[j]*sum;
01438
01439
01440
01441 if( (RateInter >= save.punarg[ipPun][1]) && (flxcor > SMALLFLOAT) )
01442 {
01443 lgPrt = true;
01444
01445 fprintf( save.ipPnunit[ipPun],
01446 "%li\t%.5e\t%.2e\t%.2e\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2e\t%.2e\t%.4s\t%.4s\n",
01447 j,
01448 AnuUnit(rfield.AnuOrg[j]),
01449 flxcor,
01450 flxcor*opac.opacity_abs[j],
01451 rfield.flux[0][j]/flxcor,
01452 rfield.otslin[j]/flxcor,
01453 rfield.otscon[j]/flxcor,
01454 (rfield.outlin[0][j] + rfield.outlin_noplot[j])/flxcor,
01455 rfield.ConInterOut[j]/flxcor,
01456 RateInter,
01457 fsum*sum,
01458 rfield.chLineLabel[j],
01459 rfield.chContLabel[j] );
01460 }
01461 }
01462 if( !lgPrt )
01463 {
01464
01465 fprintf(save.ipPnunit[ipPun],
01466 " punchdo, the PUNCH IONIZING CONTINUUM command "
01467 "did not find a strong point, sum and fsum were %.2e %.2e\n",
01468 sum,fsum);
01469 fprintf(save.ipPnunit[ipPun],
01470 " punchdo, the low-frequency energy was %.5e Ryd\n",
01471 rfield.anu[i1-1]);
01472 fprintf(save.ipPnunit[ipPun],
01473 " You can reset the threshold for the lowest fractional "
01474 "interaction to print with the second number of the save command\n"
01475 " The fraction was %.3f and this was too large.\n",
01476 save.punarg[ipPun][1]);
01477 }
01478 }
01479 }
01480 #ifdef USE_NLTE7
01481 else if( strcmp(save.chSave[ipPun],"CONl") == 0 )
01482 {
01483 if( strcmp(chTime,"LAST") != 0 )
01484 {
01485 int nEmType = (int)save.punarg[ipPun][0];
01486 ASSERT( nEmType==0 || nEmType==1 );
01487
01488 for( j=0; j<rfield.nflux; j = j+save.ncSaveSkip)
01489 {
01490
01491
01492 realnum resolution = ( save.Resolution < realnum(0.) ) ?
01493 rfield.anu[j]/rfield.widflx[j] : save.Resolution;
01494
01495
01496 flxref = (rfield.anu2[j]*((double)rfield.ConRefIncid[nEmType][j]+rfield.ConEmitReflec[nEmType][j])/rfield.widflx[j] +
01497 rfield.anu[j]*resolution*rfield.reflin[nEmType][j])*EN1RYD;
01498
01499
01500 flxatt = rfield.flux[nEmType][j]*rfield.anu2[j]*EN1RYD/
01501 rfield.widflx[j]*radius.r1r0sq *
01502 PrettyTransmission( j, rfield.trans_coef_total[j] );
01503
01504
01505 conem = (((double)rfield.ConEmitOut[nEmType][j])/
01506 rfield.widflx[j]*rfield.anu2[j] + resolution*
01507 rfield.outlin[nEmType][j]*rfield.anu[j])*radius.r1r0sq*
01508 EN1RYD*geometry.covgeo;
01509
01510
01511 flxtrn = conem + flxatt;
01512
01513
01514 double lowlim, highlim, NRGeV;
01515 lowlim = 1.0;
01516 highlim = 250.0;
01517 NRGeV = AnuUnit(rfield.AnuOrg[j]);
01518
01519 if( NRGeV >= lowlim && NRGeV <= highlim )
01520 {
01521
01522 fprintf( save.ipPnunit[ipPun],"%14.8e ", NRGeV );
01523
01524 fprintf( save.ipPnunit[ipPun],"%14.8e %14.8e %14.8e ",0.,0.,0.);
01525
01526 fprintf( save.ipPnunit[ipPun],"%14.8e\n", (flxref + flxtrn)/NRGeV );
01527 }
01528 }
01529 }
01530
01531 }
01532 #endif
01533
01534
01535 else if( strcmp(save.chSave[ipPun],"CONS") == 0 )
01536 {
01537 if( strcmp(chTime,"LAST") != 0 )
01538 {
01539
01540
01541 if( save.ipEmisFreq < 0 )
01542 save.ipEmisFreq = ipoint(save.emisfreq.Ryd());
01543 j = save.ipEmisFreq-1;
01544
01545 fprintf( save.ipPnunit[ipPun],
01546 "%.14e\t%.14e\t%.5e\t%.5e\t%.5e\n",
01547 radius.Radius_mid_zone,
01548 radius.depth_mid_zone,
01549 rfield.anu2[j]*rfield.ConEmitLocal[nzone][j]/rfield.widflx[j]*EN1RYD,
01550 opac.opacity_abs[j],
01551 opac.opacity_sct[j] );
01552 }
01553 }
01554
01555 else if( strcmp(save.chSave[ipPun],"CORA") == 0 )
01556 {
01557
01558
01559
01560
01561 if( strcmp(chTime,"LAST") == 0 )
01562 {
01563
01564 for( j=0;j<rfield.nflux;j = j + save.ncSaveSkip)
01565 {
01566 fprintf( save.ipPnunit[ipPun],
01567 "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%4.4s\t%4.4s\t",
01568 AnuUnit(rfield.AnuOrg[j]),
01569 rfield.flux[0][j],
01570 rfield.otslin[j],
01571 rfield.otscon[j],
01572 rfield.ConRefIncid[0][j],
01573 rfield.ConEmitReflec[0][j],
01574 rfield.ConInterOut[j],
01575 rfield.outlin[0][j]+rfield.outlin_noplot[j],
01576 rfield.ConEmitOut[0][j],
01577 rfield.chLineLabel[j],
01578 rfield.chContLabel[j]
01579 );
01580
01581 fprintf( save.ipPnunit[ipPun], "%li\n", rfield.line_count[j] );
01582 }
01583 }
01584 }
01585
01586 else if( strcmp(save.chSave[ipPun],"CONo") == 0 )
01587 {
01588
01589
01590
01591
01592 if( strcmp(chTime,"LAST") == 0 )
01593 {
01594
01595 for( j=0;j<rfield.nflux; j = j + save.ncSaveSkip)
01596 {
01597 fprintf( save.ipPnunit[ipPun], "%11.5e%10.2e%10.2e\n",
01598 AnuUnit(rfield.AnuOrg[j]),
01599 rfield.ConEmitOut[0][j]+ rfield.outlin[0][j] + rfield.outlin_noplot[j],
01600 (rfield.flux[0][j] + rfield.otscon[j] + rfield.otslin[j] +
01601 rfield.ConInterOut[j])*opac.opacity_abs[j]*
01602 rfield.anu[j] );
01603 }
01604 }
01605 }
01606
01607 else if( strcmp(save.chSave[ipPun],"CONO") == 0 )
01608 {
01609
01610
01611
01612
01613 if( strcmp(chTime,"LAST") == 0 )
01614 {
01615
01616 for( j=0; j<rfield.nflux; j = j + save.ncSaveSkip)
01617 {
01618
01619 realnum resolution = ( save.Resolution < realnum(0.) ) ?
01620 rfield.anu[j]/rfield.widflx[j] : save.Resolution;
01621
01622 fprintf( save.ipPnunit[ipPun], "%11.5e%10.2e%10.2e%10.2e%10.2e\n",
01623 AnuUnit(rfield.AnuOrg[j]),
01624 rfield.flux[0][j]*rfield.anu2[j]* EN1RYD/rfield.widflx[j]*radius.r1r0sq,
01625 rfield.ConInterOut[j]/rfield.widflx[j]*rfield.anu2[j]*radius.r1r0sq*EN1RYD,
01626 resolution*(rfield.outlin[0][j]+rfield.outlin_noplot[j])*rfield.anu[j]*radius.r1r0sq*EN1RYD,
01627 (rfield.ConInterOut[j]/rfield.widflx[j]*
01628 rfield.anu2[j] + resolution*(rfield.outlin[0][j]+rfield.outlin_noplot[j])*
01629 rfield.anu[j])*radius.r1r0sq*EN1RYD );
01630 }
01631 }
01632 }
01633
01634 else if( strcmp(save.chSave[ipPun],"CONT") == 0 )
01635 {
01636
01637
01638
01639
01640
01641 if( strcmp(chTime,"LAST") == 0 )
01642 {
01643 fprintf( save.ipPnunit[ipPun], "#\n" );
01644 fprintf( save.ipPnunit[ipPun], "%32ld # file format version number\n",
01645 VERSION_TRNCON );
01646 fprintf( save.ipPnunit[ipPun], "%s # check 1\n",
01647 continuum.mesh_md5sum.c_str() );
01648 union {
01649 double x;
01650 uint32 i[2];
01651 } u;
01652 u.x = double(rfield.emm);
01653 if( cpu.i().big_endian() )
01654 fprintf( save.ipPnunit[ipPun], "%23.8x %8.8x # check 2\n",
01655 u.i[0], u.i[1] );
01656 else
01657 fprintf( save.ipPnunit[ipPun], "%23.8x %8.8x # check 2\n",
01658 u.i[1], u.i[0] );
01659 u.x = double(rfield.egamry);
01660 if( cpu.i().big_endian() )
01661 fprintf( save.ipPnunit[ipPun], "%23.8x %8.8x # check 3\n",
01662 u.i[0], u.i[1] );
01663 else
01664 fprintf( save.ipPnunit[ipPun], "%23.8x %8.8x # check 3\n",
01665 u.i[1], u.i[0] );
01666 u.x = continuum.ResolutionScaleFactor;
01667 if( cpu.i().big_endian() )
01668 fprintf( save.ipPnunit[ipPun], "%23.8x %8.8x # check 4\n",
01669 u.i[0], u.i[1] );
01670 else
01671 fprintf( save.ipPnunit[ipPun], "%23.8x %8.8x # check 4\n",
01672 u.i[1], u.i[0] );
01673 fprintf( save.ipPnunit[ipPun], "%32ld # nflux\n",
01674 (rfield.nflux+save.ncSaveSkip-1)/save.ncSaveSkip );
01675 fprintf( save.ipPnunit[ipPun], "#\n" );
01676
01677 const realnum *trans_coef_total=rfield.getCoarseTransCoef();
01678
01679
01680 for( j=0; j < rfield.nflux; j += save.ncSaveSkip )
01681 {
01682
01683
01684
01685
01686 flxatt = rfield.flux[0][j]*rfield.anu2[j]*EN1RYD/
01687 rfield.widflx[j]*radius.r1r0sq*trans_coef_total[j];
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704 conem = (rfield.ConEmitOut[0][j] + rfield.outlin[0][j]) / rfield.widflx[j]*
01705 rfield.anu2[j]*radius.r1r0sq*EN1RYD*geometry.covgeo;
01706
01707 flxtrn = conem + flxatt;
01708
01709
01710
01711
01712 fprintf( save.ipPnunit[ipPun],"%.5e\t%.3e\t%.3e\n",
01713 AnuUnit(rfield.AnuOrg[j]), flxtrn,
01714 trans_coef_total[j] );
01715 }
01716 }
01717 }
01718
01719 else if( strcmp(save.chSave[ipPun],"CON2") == 0 )
01720 {
01721
01722 if( strcmp(chTime,"LAST") == 0 )
01723 {
01724
01725 for( j=0; j<rfield.nflux; j = j+save.ncSaveSkip)
01726 {
01727 fprintf( save.ipPnunit[ipPun], "%.5e\t%.5e\t%.5e\n",
01728 AnuUnit(rfield.AnuOrg[j]),
01729 rfield.TotDiff2Pht[j]/rfield.widflx[j] ,
01730 rfield.TotDiff2Pht[j]*rfield.anu2[j]*EN1RYD/rfield.widflx[j]);
01731 }
01732 }
01733 }
01734
01735 else if( strcmp(save.chSave[ipPun],"DUSE") == 0 )
01736 {
01737
01738 if( strcmp(chTime,"LAST") != 0 )
01739 {
01740 fprintf( save.ipPnunit[ipPun], " %.5e\t",
01741 radius.depth_mid_zone );
01742
01743
01744 fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended);
01745
01746
01747 fprintf( save.ipPnunit[ipPun], "%.2e\n" , rfield.extin_mag_V_point);
01748 }
01749 }
01750
01751 else if( strcmp(save.chSave[ipPun],"DUSO") == 0 )
01752 {
01753
01754 if( strcmp(chTime,"LAST") == 0 )
01755 {
01756 for( j=0; j < rfield.nflux; j++ )
01757 {
01758 double scat;
01759 fprintf( save.ipPnunit[ipPun],
01760 "%.5e\t%.2e\t%.2e\t%.2e\t",
01761
01762 AnuUnit(rfield.AnuOrg[j]),
01763
01764 gv.dstab[j] + gv.dstsc[j],
01765
01766 gv.dstab[j],
01767
01768 gv.dstsc[j] );
01769
01770 scat = 0.;
01771
01772 for( size_t nd=0; nd < gv.bin.size(); nd++ )
01773 {
01774 scat += gv.bin[nd]->pure_sc1[j]*gv.bin[nd]->dstAbund;
01775 }
01776
01777 fprintf( save.ipPnunit[ipPun],
01778 "%.2e\t", scat );
01779 fprintf( save.ipPnunit[ipPun],
01780 "%.2e\n", gv.dstsc[j]/SDIV(gv.dstab[j] + gv.dstsc[j]) );
01781 }
01782 }
01783 }
01784
01785
01786 else if( strcmp(save.chSave[ipPun],"DUSA") == 0 ||
01787 strcmp(save.chSave[ipPun],"DUSD") == 0 )
01788 {
01789 bool lgDGRatio = ( strcmp(save.chSave[ipPun],"DUSD") == 0 );
01790
01791
01792 if( strcmp(chTime,"LAST") != 0 )
01793 {
01794
01795 static bool lgMustPrtHeaderDRRatio = true,
01796 lgMustPrtHeaderAbundance=true;
01797
01798 if( ( lgMustPrtHeaderDRRatio && lgDGRatio ) ||
01799 ( lgMustPrtHeaderAbundance && !lgDGRatio ) )
01800 {
01801
01802
01803 if( lgMustPrtHeaderDRRatio && lgDGRatio )
01804 lgMustPrtHeaderDRRatio = false;
01805 else if( lgMustPrtHeaderAbundance &&!lgDGRatio )
01806 lgMustPrtHeaderAbundance = false;
01807
01808 fprintf( save.ipPnunit[ipPun], "#Depth" );
01809 for( size_t nd=0; nd < gv.bin.size(); ++nd )
01810 fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
01811 fprintf( save.ipPnunit[ipPun], "\ttotal\n" );
01812
01813
01814 fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" );
01815 for( size_t nd=0; nd < gv.bin.size(); ++nd )
01816 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
01817 fprintf( save.ipPnunit[ipPun], "\txxxx\n" );
01818 }
01819 fprintf( save.ipPnunit[ipPun], " %.5e",
01820 radius.depth_mid_zone );
01821
01822 double total = 0.;
01823 for( size_t nd=0; nd < gv.bin.size(); ++nd )
01824 {
01825 double abund = gv.bin[nd]->IntVol*gv.bin[nd]->dustp[0]*
01826 gv.bin[nd]->cnv_H_pCM3;
01827 if( lgDGRatio )
01828 abund /= dense.xMassDensity;
01829 fprintf( save.ipPnunit[ipPun], "\t%.3e", abund );
01830 total += abund;
01831 }
01832 fprintf( save.ipPnunit[ipPun], "\t%.3e\n", total );
01833 }
01834 }
01835
01836 else if( strcmp(save.chSave[ipPun],"DUSP") == 0 )
01837 {
01838
01839 if( strcmp(chTime,"LAST") != 0 )
01840 {
01841
01842 static bool lgMustPrtHeader = true;
01843
01844 if( lgMustPrtHeader )
01845 {
01846
01847 fprintf( save.ipPnunit[ipPun], "#Depth" );
01848 for( size_t nd=0; nd < gv.bin.size(); ++nd )
01849 fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
01850 fprintf( save.ipPnunit[ipPun], "\n" );
01851
01852
01853 fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" );
01854 for( size_t nd=0; nd < gv.bin.size(); ++nd )
01855 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
01856 fprintf( save.ipPnunit[ipPun], "\n" );
01857
01858
01859 lgMustPrtHeader = false;
01860 }
01861 fprintf( save.ipPnunit[ipPun], " %.5e",
01862 radius.depth_mid_zone );
01863
01864 for( size_t nd=0; nd < gv.bin.size(); ++nd )
01865 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->dstpot*EVRYD );
01866 fprintf( save.ipPnunit[ipPun], "\n" );
01867 }
01868 }
01869
01870 else if( strcmp(save.chSave[ipPun],"DUSR") == 0 )
01871 {
01872
01873 if( strcmp(chTime,"LAST") != 0 )
01874 {
01875
01876 static bool lgMustPrtHeader = true;
01877
01878 if( lgMustPrtHeader )
01879 {
01880
01881 fprintf( save.ipPnunit[ipPun], "#Depth" );
01882 for( size_t nd=0; nd < gv.bin.size(); ++nd )
01883 fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
01884 fprintf( save.ipPnunit[ipPun], "\n" );
01885
01886
01887 fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" );
01888 for( size_t nd=0; nd < gv.bin.size(); ++nd )
01889 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
01890 fprintf( save.ipPnunit[ipPun], "\n" );
01891
01892
01893 lgMustPrtHeader = false;
01894 }
01895 fprintf( save.ipPnunit[ipPun], " %.5e",
01896 radius.depth_mid_zone );
01897
01898 for( size_t nd=0; nd < gv.bin.size(); ++nd )
01899 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->rate_h2_form_grains_used );
01900 fprintf( save.ipPnunit[ipPun], "\n" );
01901 }
01902 }
01903
01904 else if( strcmp(save.chSave[ipPun],"DUST") == 0 )
01905 {
01906
01907 if( strcmp(chTime,"LAST") != 0 )
01908 {
01909
01910 static bool lgMustPrtHeader = true;
01911
01912 if( lgMustPrtHeader )
01913 {
01914
01915 fprintf( save.ipPnunit[ipPun], "#Depth" );
01916 for( size_t nd=0; nd < gv.bin.size(); ++nd )
01917 fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
01918 fprintf( save.ipPnunit[ipPun], "\n" );
01919
01920
01921 fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" );
01922 for( size_t nd=0; nd < gv.bin.size(); ++nd )
01923 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
01924 fprintf( save.ipPnunit[ipPun], "\n" );
01925
01926
01927 lgMustPrtHeader = false;
01928 }
01929 fprintf( save.ipPnunit[ipPun], " %.5e",
01930 radius.depth_mid_zone );
01931 for( size_t nd=0; nd < gv.bin.size(); ++nd )
01932 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->tedust );
01933 fprintf( save.ipPnunit[ipPun], "\n" );
01934 }
01935 }
01936
01937 else if( strcmp(save.chSave[ipPun],"DUSC") == 0 )
01938 {
01939
01940
01941 if( strcmp(chTime,"LAST") != 0 )
01942 {
01943
01944 static bool lgMustPrtHeader = true;
01945
01946 if( lgMustPrtHeader )
01947 {
01948
01949 fprintf( save.ipPnunit[ipPun], "#Depth\tne(grn)" );
01950 for( size_t nd=0; nd < gv.bin.size(); ++nd )
01951 fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
01952 fprintf( save.ipPnunit[ipPun], "\n" );
01953
01954
01955 fprintf( save.ipPnunit[ipPun], "#grn rad (mic)\tne\t" );
01956 for( size_t nd=0; nd < gv.bin.size(); ++nd )
01957 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
01958 fprintf( save.ipPnunit[ipPun], "\n" );
01959
01960
01961 lgMustPrtHeader = false;
01962 }
01963
01964 fprintf( save.ipPnunit[ipPun], " %.5e\t%.4e",
01965 radius.depth_mid_zone ,
01966
01967
01968 gv.TotalEden );
01969
01970
01971 for( size_t nd=0; nd < gv.bin.size(); ++nd )
01972 {
01973 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AveDustZ );
01974 }
01975 fprintf( save.ipPnunit[ipPun], "\n" );
01976 }
01977 }
01978
01979 else if( strcmp(save.chSave[ipPun],"DUSH") == 0 )
01980 {
01981
01982 if( strcmp(chTime,"LAST") != 0 )
01983 {
01984
01985 static bool lgMustPrtHeader = true;
01986
01987 if( lgMustPrtHeader )
01988 {
01989
01990 fprintf( save.ipPnunit[ipPun], "#Depth" );
01991 for( size_t nd=0; nd < gv.bin.size(); ++nd )
01992 fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
01993 fprintf( save.ipPnunit[ipPun], "\n" );
01994
01995
01996 fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" );
01997 for( size_t nd=0; nd < gv.bin.size(); ++nd )
01998 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
01999 fprintf( save.ipPnunit[ipPun], "\n" );
02000
02001
02002 lgMustPrtHeader = false;
02003 }
02004 fprintf( save.ipPnunit[ipPun], " %.5e",
02005 radius.depth_mid_zone );
02006
02007 for( size_t nd=0; nd < gv.bin.size(); ++nd )
02008 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->GasHeatPhotoEl );
02009 fprintf( save.ipPnunit[ipPun], "\n" );
02010 }
02011 }
02012
02013 else if( strcmp(save.chSave[ipPun],"DUSV") == 0 )
02014 {
02015
02016 if( strcmp(chTime,"LAST") != 0 )
02017 {
02018
02019 static bool lgMustPrtHeader = true;
02020
02021 if( lgMustPrtHeader )
02022 {
02023
02024 fprintf( save.ipPnunit[ipPun], "#Depth" );
02025 for( size_t nd=0; nd < gv.bin.size(); ++nd )
02026 fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
02027 fprintf( save.ipPnunit[ipPun], "\n" );
02028
02029
02030 fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" );
02031 for( size_t nd=0; nd < gv.bin.size(); ++nd )
02032 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
02033 fprintf( save.ipPnunit[ipPun], "\n" );
02034
02035
02036 lgMustPrtHeader = false;
02037 }
02038 fprintf( save.ipPnunit[ipPun], " %.5e",
02039 radius.depth_mid_zone );
02040
02041 for( size_t nd=0; nd < gv.bin.size(); ++nd )
02042 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->DustDftVel*1e-5 );
02043 fprintf( save.ipPnunit[ipPun], "\n" );
02044 }
02045 }
02046
02047
02048 else if( strcmp(save.chSave[ipPun],"DUSQ") == 0 )
02049 {
02050
02051 if( strcmp(chTime,"LAST") == 0 )
02052 {
02053 for( j=0; j < rfield.nflux; j++ )
02054 {
02055 fprintf( save.ipPnunit[ipPun], "%11.4e",
02056 rfield.anu[j] );
02057 for( size_t nd=0; nd < gv.bin.size(); nd++ )
02058 {
02059 fprintf( save.ipPnunit[ipPun], "%9.1e%9.1e",
02060 gv.bin[nd]->dstab1[j]*4./gv.bin[nd]->IntArea,
02061 gv.bin[nd]->pure_sc1[j]*gv.bin[nd]->asym[j]*4./gv.bin[nd]->IntArea );
02062 }
02063 fprintf( save.ipPnunit[ipPun], "\n" );
02064 }
02065 }
02066 }
02067
02068 else if( strcmp(save.chSave[ipPun],"ELEM") == 0 )
02069 {
02070 if( strcmp(chTime,"LAST") != 0 )
02071 {
02072 realnum renorm = 1.f;
02073
02074
02075
02076 nelem = (long int)save.punarg[ipPun][0];
02077 ASSERT( nelem >= ipHYDROGEN );
02078
02079
02080 if( dense.lgElmtOn[nelem] )
02081 {
02082
02083
02084 if( save.punarg[ipPun][1] == 0 )
02085 renorm = dense.gas_phase[nelem];
02086
02087 fprintf( save.ipPnunit[ipPun], " %.5e", radius.depth_mid_zone );
02088
02089 for( j=0; j <= (nelem + 1); ++j)
02090 {
02091 fprintf( save.ipPnunit[ipPun], "\t%.2e",
02092 dense.xIonDense[nelem][j]/renorm );
02093 }
02094 if( nelem==ipHYDROGEN )
02095 {
02096
02097 fprintf( save.ipPnunit[ipPun], "\t%.2e",
02098 hmi.H2_total/renorm );
02099 }
02100
02101 else if( nelem==ipCARBON )
02102 {
02103 fprintf( save.ipPnunit[ipPun], "\t%.2e\t%.2e\t%.2e",
02104 colden.C1Pops[0]/renorm, colden.C1Pops[1]/renorm, colden.C1Pops[2]/renorm);
02105 fprintf( save.ipPnunit[ipPun], "\t%.2e\t%.2e",
02106 colden.C2Pops[0]/renorm, colden.C2Pops[1]/renorm);
02107 fprintf( save.ipPnunit[ipPun], "\t%.2e",
02108 findspecieslocal("CO")->den/renorm );
02109 }
02110 else if( nelem==ipOXYGEN )
02111 {
02112 fprintf( save.ipPnunit[ipPun], "\t%.2e\t%.2e\t%.2e",
02113 colden.O1Pops[0]/renorm, colden.O1Pops[1]/renorm, colden.O1Pops[2]/renorm);
02114 }
02115 fprintf( save.ipPnunit[ipPun], "\n" );
02116 }
02117 }
02118 }
02119
02120 else if( strcmp(save.chSave[ipPun],"RECA") == 0 )
02121 {
02122
02123
02124 ion_recombAGN( save.ipPnunit[ipPun] );
02125 cdEXIT(EXIT_FAILURE);
02126 }
02127
02128 else if( strcmp(save.chSave[ipPun],"RECE") == 0 )
02129 {
02130
02131
02132
02133
02134
02135
02136 fprintf( save.ipPnunit[ipPun],
02137 "%12.4e %12.4e %12.4e %12.4e\n",
02138 iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].RadRecomb[ipRecRad],
02139 iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].RadRecomb[ipRecNetEsc] ,
02140 iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].RadRecomb[ipRecRad],
02141 iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].RadRecomb[ipRecNetEsc]);
02142 }
02143
02144 else if( strncmp(save.chSave[ipPun],"FEc" , 3 ) == 0 )
02145 {
02146
02147 if( strcmp(chTime,"LAST") == 0 )
02148 {
02149
02150 if( save.punarg[ipPun][0] == 1 )
02151 {
02152
02153 char chTempHeader[100];
02154 long ipFeII_Cont_type;
02155 if( strcmp(save.chSave[ipPun],"FEcI") == 0 )
02156 {
02157 strcpy(chTempHeader , "#FeII inward: Wl(A)\tInt[erg cm-2 s-1]\n");
02158 ipFeII_Cont_type = 1;
02159 }
02160 else if( strcmp(save.chSave[ipPun],"FEcO") == 0 )
02161 {
02162 strcpy(chTempHeader , "#FeII outward: Wl(A)\tInt[erg cm-2 s-1]\n");
02163 ipFeII_Cont_type = 2;
02164 }
02165 else if( strcmp(save.chSave[ipPun],"FEcT") == 0 )
02166 {
02167 strcpy(chTempHeader , "#FeII total: Wl(A)\tInt[erg cm-2 s-1]\n");
02168 ipFeII_Cont_type = 3;
02169 }
02170 else
02171 TotalInsanity();
02172
02173 if( save.lgPunHeader[ipPun] &&
02174 (!grid.lgGrid || optimize.nOptimiz==0) )
02175 {
02176
02177 j = 0;
02178 fprintf( save.ipPnunit[ipPun], "%s%.5f",
02179 chTempHeader ,
02180 AnuUnit((realnum)RYDLAM/((FeII_Cont[j][0]+FeII_Cont[j+1][0])/2.f) ));
02181 for( j=1; j < nFeIIConBins; j++ )
02182 {
02183 fprintf( save.ipPnunit[ipPun], "\t%.5e",
02184 AnuUnit((realnum)RYDLAM/((FeII_Cont[j][0]+FeII_Cont[j+1][0])/2.f) ));
02185 }
02186 fprintf( save.ipPnunit[ipPun], "\n");
02187 save.lgPunHeader[ipPun] = false;
02188 }
02189
02190 fprintf( save.ipPnunit[ipPun], "%.2f",
02191 SaveFeII_cont( 0 , ipFeII_Cont_type ));
02192 for( j=1; j < nFeIIConBins; j++ )
02193 {
02194 fprintf( save.ipPnunit[ipPun], "\t%e",
02195 SaveFeII_cont( j , ipFeII_Cont_type ) );
02196 }
02197 fprintf( save.ipPnunit[ipPun], "\n");
02198 }
02199 else if( save.punarg[ipPun][0] == 2 )
02200 {
02201
02202
02203 if( save.lgPunHeader[ipPun] &&
02204 (!grid.lgGrid || optimize.nOptimiz==0) )
02205 {
02206 fprintf( save.ipPnunit[ipPun], "FeiI wl, total, inward, outward\n");
02207 save.lgPunHeader[ipPun] = false;
02208 }
02209 for( j=0; j < nFeIIConBins; j++ )
02210 {
02211 fprintf( save.ipPnunit[ipPun], "%.5f",
02212 AnuUnit((realnum)RYDLAM/((FeII_Cont[j][0]+FeII_Cont[j+1][0])/2.f) ));
02213 fprintf( save.ipPnunit[ipPun], "\t%e",
02214 SaveFeII_cont( j , 3 ) );
02215 fprintf( save.ipPnunit[ipPun], "\t%e",
02216 SaveFeII_cont( j , 1 ) );
02217 fprintf( save.ipPnunit[ipPun], "\t%e",
02218 SaveFeII_cont( j , 2 ) );
02219 fprintf( save.ipPnunit[ipPun], "\n");
02220 }
02221 }
02222 }
02223 }
02224
02225
02226 else if( strcmp(save.chSave[ipPun],"FENl") == 0 )
02227 {
02228 if( strcmp(chTime,"LAST") == 0 )
02229 {
02230
02231 FeIIPunchColden( save.ipPnunit[ipPun] );
02232 }
02233 }
02234
02235 else if( strcmp(save.chSave[ipPun],"FE2l") == 0 )
02236 {
02237 if( strcmp(chTime,"LAST") == 0 )
02238 {
02239
02240 FeIIPunchLevels( save.ipPnunit[ipPun] );
02241 }
02242 }
02243
02244 else if( strcmp(save.chSave[ipPun],"FEli") == 0 )
02245 {
02246 if( strcmp(chTime,"LAST") == 0 )
02247 {
02248
02249 FeIISaveLines( save.ipPnunit[ipPun] );
02250 }
02251 }
02252
02253 else if( strcmp(save.chSave[ipPun],"FE2o") == 0 )
02254 {
02255 if( strcmp(chTime,"LAST") == 0 )
02256 {
02257
02258 FeIIPunchOpticalDepth( save.ipPnunit[ipPun] );
02259 }
02260 }
02261
02262 else if( strcmp(save.chSave[ipPun],"FRED") == 0 )
02263 {
02264
02265
02266 if( strcmp(chTime,"LAST") != 0 )
02267 {
02268
02269 fprintf( save.ipPnunit[ipPun], "%.5e\t%.5e\t%.3e\t%.3e\t%.3e"
02270 "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e"
02271 "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e"
02272 "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e"
02273 "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e"
02274 "\t%.3e\t%.3e\n",
02275 radius.Radius, radius.depth ,wind.windv/1e5,
02276 wind.dvdr,
02277 dense.gas_phase[ipHYDROGEN], dense.eden , phycon.te,
02278 wind.AccelLine , wind.AccelCont ,
02279 wind.fmul ,
02280
02281
02282 pressure.pinzon_PresIntegElecThin/dense.xMassDensity/radius.drad_x_fillfac ,
02283 mean.xIonMean[0][ipHYDROGEN][0][0] , mean.xIonMean[0][ipHYDROGEN][1][0] ,
02284 mean.xIonMean[0][ipHELIUM][0][0] , mean.xIonMean[0][ipHELIUM][1][0] ,
02285 mean.xIonMean[0][ipHELIUM][2][0] ,
02286 mean.xIonMean[0][ipCARBON][1][0] , mean.xIonMean[0][ipCARBON][2][0] ,
02287 mean.xIonMean[0][ipCARBON][3][0] ,
02288 mean.xIonMean[0][ipOXYGEN][0][0] , mean.xIonMean[0][ipOXYGEN][1][0] ,
02289 mean.xIonMean[0][ipOXYGEN][2][0] , mean.xIonMean[0][ipOXYGEN][3][0] ,
02290 mean.xIonMean[0][ipOXYGEN][4][0] , mean.xIonMean[0][ipOXYGEN][5][0] ,
02291 mean.xIonMean[0][ipOXYGEN][6][0] , mean.xIonMean[0][ipOXYGEN][7][0] ,
02292 dense.xIonDense[ipHYDROGEN][0] , dense.xIonDense[ipHYDROGEN][1] ,
02293 dense.xIonDense[ipHELIUM][0] , dense.xIonDense[ipHELIUM][1] ,
02294 dense.xIonDense[ipHELIUM][2] ,
02295 dense.xIonDense[ipCARBON][1] , dense.xIonDense[ipCARBON][2] ,
02296 dense.xIonDense[ipCARBON][3] ,
02297 dense.xIonDense[ipOXYGEN][0] , dense.xIonDense[ipOXYGEN][1] ,
02298 dense.xIonDense[ipOXYGEN][2] , dense.xIonDense[ipOXYGEN][3] ,
02299 dense.xIonDense[ipOXYGEN][4] , dense.xIonDense[ipOXYGEN][5] ,
02300 dense.xIonDense[ipOXYGEN][6] , dense.xIonDense[ipOXYGEN][7] ,
02301 mean.xIonMean[0][ipMAGNESIUM][1][0] , dense.xIonDense[ipMAGNESIUM][1],
02302 TauLines[ipT1032].Emis().TauIn() ,
02303 TauLines[ipT1032].Emis().TauCon() );
02304 }
02305 }
02306
02307 else if( strcmp(save.chSave[ipPun],"FE2d") == 0 )
02308 {
02309
02310 if( strcmp(chTime,"LAST") != 0 )
02311 FeIIPunDepart(save.ipPnunit[ipPun],false);
02312 }
02313
02314 else if( strcmp(save.chSave[ipPun],"FE2D") == 0 )
02315 {
02316
02317 if( strcmp(chTime,"LAST") != 0 )
02318 FeIIPunDepart(save.ipPnunit[ipPun],true);
02319 }
02320
02321 else if( strcmp(save.chSave[ipPun],"FE2p") == 0 )
02322 {
02323 bool lgFlag = false;
02324 if( save.punarg[ipPun][2] )
02325 lgFlag = true;
02326
02327 if( strcmp(chTime,"LAST") != 0 )
02328 FeIIPunPop(save.ipPnunit[ipPun],false,0,0,
02329 lgFlag);
02330 }
02331
02332 else if( strcmp(save.chSave[ipPun],"FE2P") == 0 )
02333 {
02334 bool lgFlag = false;
02335 if( save.punarg[ipPun][2] )
02336 lgFlag = true;
02337
02338 if( strcmp(chTime,"LAST") != 0 )
02339 FeIIPunPop(save.ipPnunit[ipPun],
02340 true,
02341 (long int)save.punarg[ipPun][0],
02342 (long int)save.punarg[ipPun][1],
02343 lgFlag);
02344 }
02345
02346
02347 else if( strcmp(save.chSave[ipPun],"FITS") == 0 )
02348 {
02349 if( strcmp(chTime,"LAST") == 0 )
02350 {
02351 saveFITSfile( save.ipPnunit[ipPun], NUM_OUTPUT_TYPES );
02352 }
02353 }
02354
02355 else if( strcmp(save.chSave[ipPun],"GAMt") == 0 )
02356 {
02357 if( strcmp(chTime,"LAST") != 0 )
02358 {
02359 long ns;
02360
02361 for( nelem=0; nelem < LIMELM; nelem++ )
02362 {
02363 for( ion=0; ion <= nelem; ion++ )
02364 {
02365 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
02366 {
02367 fprintf( save.ipPnunit[ipPun], "%3ld%3ld%3ld%10.2e%10.2e%10.2e",
02368 nelem+1, ion+1, ns+1,
02369 ionbal.PhotoRate_Shell[nelem][ion][ns][0],
02370 ionbal.PhotoRate_Shell[nelem][ion][ns][1] ,
02371 ionbal.PhotoRate_Shell[nelem][ion][ns][2] );
02372
02373 for( j=0; j < t_yield::Inst().nelec_eject(nelem,ion,ns); j++ )
02374 {
02375 fprintf( save.ipPnunit[ipPun], "%5.2f",
02376 t_yield::Inst().elec_eject_frac(nelem,ion,ns,j) );
02377 }
02378 fprintf( save.ipPnunit[ipPun], "\n" );
02379 }
02380 }
02381 }
02382 }
02383 }
02384
02385
02386 else if( strcmp(save.chSave[ipPun],"GAMe") == 0 )
02387 {
02388 if( strcmp(chTime,"LAST") != 0 )
02389 {
02390 int ns;
02391 nelem = (long)save.punarg[ipPun][0];
02392 ion = (long)save.punarg[ipPun][1];
02393
02394 ns = Heavy.nsShells[nelem][ion]-1;
02395
02396 GammaPrt(
02397 opac.ipElement[nelem][ion][ns][0] ,
02398 opac.ipElement[nelem][ion][ns][1] ,
02399 opac.ipElement[nelem][ion][ns][2] ,
02400 save.ipPnunit[ipPun],
02401 ionbal.PhotoRate_Shell[nelem][ion][ns][0] ,
02402 ionbal.PhotoRate_Shell[nelem][ion][ns][0]*0.1 );
02403 }
02404 }
02405
02406 else if( strcmp(save.chSave[ipPun],"GAUN") == 0 )
02407 {
02408
02409 if( strcmp(chTime,"LAST") != 0 )
02410 SaveGaunts(save.ipPnunit[ipPun]);
02411 }
02412
02413
02414
02415 }
02416 else
02417 {
02418
02419 lgNoHitFirstBranch = true;
02420 }
02421
02422
02423
02424 if( iterations.lgLastIt || !save.lgPunLstIter[ipPun] ||
02425 ( lgAbort && strcmp(chTime,"LAST") == 0 ) ||
02426 ( dynamics.lgTimeDependentStatic && dynamics.lgStatic_completed ) )
02427 {
02428 if( strcmp(save.chSave[ipPun],"HISp") == 0 )
02429 {
02430
02431 if( strcmp(chTime,"LAST") != 0 )
02432 {
02433
02434 if( !conv.lgConvPres )
02435 {
02436 fprintf( save.ipPnunit[ipPun],
02437 "#PROBLEM Pressure not converged iter %li zone %li density-pressure follows:\n",
02438 iteration , nzone );
02439 }
02440
02441 if( !conv.lgConvTemp )
02442 {
02443 fprintf( save.ipPnunit[ipPun],
02444 "#PROBLEM Temperature not converged iter %li zone %li density-pressure follows:\n",
02445 iteration , nzone );
02446 }
02447 for( unsigned long k=0; k < conv.hist_pres_density.size(); ++k )
02448 {
02449
02450 fprintf( save.ipPnunit[ipPun] , "%2li %4li\t%.5e\t%.5e\t%.5e\n",
02451 iteration,
02452 nzone,
02453 conv.hist_pres_density[k],
02454 conv.hist_pres_current[k],
02455 conv.hist_pres_error[k]);
02456 }
02457 }
02458 }
02459
02460 else if( strcmp(save.chSave[ipPun],"HISt") == 0 )
02461 {
02462
02463 if( strcmp(chTime,"LAST") != 0 )
02464 {
02465
02466 if( !conv.lgConvPres )
02467 {
02468 fprintf( save.ipPnunit[ipPun],
02469 "#PROBLEM Pressure not converged iter %li zone %li temp heat cool follows:\n",
02470 iteration , nzone );
02471 }
02472
02473 if( !conv.lgConvTemp )
02474 {
02475 fprintf( save.ipPnunit[ipPun],
02476 "#PROBLEM Temperature not converged iter %li zone %li temp heat cool follows:\n",
02477 iteration , nzone );
02478 }
02479 for( unsigned long k=0; k < conv.hist_temp_temp.size(); ++k )
02480 {
02481
02482 fprintf( save.ipPnunit[ipPun] , "%2li %4li\t%.5e\t%.5e\t%.5e\n",
02483 iteration,
02484 nzone,
02485 conv.hist_temp_temp[k],
02486 conv.hist_temp_heat[k],
02487 conv.hist_temp_cool[k]);
02488 }
02489 }
02490 }
02491
02492 else if( strncmp(save.chSave[ipPun],"H2",2) == 0 )
02493 {
02494
02495 save.whichDiatomToPrint[ipPun]->H2_PunchDo( save.ipPnunit[ipPun] , save.chSave[ipPun] , chTime, ipPun );
02496 }
02497
02498 else if( strcmp(save.chSave[ipPun],"HEAT") == 0 )
02499 {
02500
02501 if( strcmp(chTime,"LAST") != 0 )
02502 SaveHeat(save.ipPnunit[ipPun]);
02503 }
02504
02505 else if( strncmp(save.chSave[ipPun],"HE",2) == 0 )
02506 {
02507
02508
02509 if( strcmp(save.chSave[ipPun] , "HELW") == 0 )
02510 {
02511 if( strcmp(chTime,"LAST") == 0 )
02512 {
02513
02514 fprintf( save.ipPnunit[ipPun],
02515 "Z\tElem\t2 1P->1 1S\t2 3P1->1 1S\t2 3P2->1 1S"
02516 "\t2 3S->1 1S\t2 3P2->2 3S\t2 3P1->2 3S\t2 3P0->2 3S" );
02517 fprintf( save.ipPnunit[ipPun], "\n" );
02518 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
02519 {
02520
02521 fprintf( save.ipPnunit[ipPun], "%li\t%s",
02522 nelem+1 , elementnames.chElementSym[nelem] );
02523
02524 fprintf( save.ipPnunit[ipPun], "\t" );
02525 prt_wl( save.ipPnunit[ipPun] ,
02526 iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).WLAng() );
02527 fprintf( save.ipPnunit[ipPun], "\t" );
02528 prt_wl( save.ipPnunit[ipPun] ,
02529 iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P1,ipHe1s1S).WLAng() );
02530 fprintf( save.ipPnunit[ipPun], "\t" );
02531 prt_wl( save.ipPnunit[ipPun] ,
02532 iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P2,ipHe1s1S).WLAng() );
02533 fprintf( save.ipPnunit[ipPun], "\t" );
02534 prt_wl( save.ipPnunit[ipPun] ,
02535 iso_sp[ipHE_LIKE][nelem].trans(ipHe2s3S,ipHe1s1S).WLAng() );
02536 fprintf( save.ipPnunit[ipPun], "\t" );
02537 prt_wl( save.ipPnunit[ipPun] ,
02538 iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P2,ipHe2s3S).WLAng() );
02539 fprintf( save.ipPnunit[ipPun], "\t" );
02540 prt_wl( save.ipPnunit[ipPun] ,
02541 iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P1,ipHe2s3S).WLAng() );
02542 fprintf( save.ipPnunit[ipPun], "\t" );
02543 prt_wl( save.ipPnunit[ipPun] ,
02544 iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P0,ipHe2s3S).WLAng() );
02545 fprintf( save.ipPnunit[ipPun], "\n");
02546 }
02547 }
02548 }
02549 else
02550 TotalInsanity();
02551 }
02552
02553
02554 else if( strcmp(save.chSave[ipPun],"HUMM") == 0 )
02555 {
02556 eps = iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Aul()/
02557 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Coll().ColUL( colliders );
02558 fprintf( save.ipPnunit[ipPun],
02559 " %.5e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",
02560 radius.depth_mid_zone,
02561 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauIn(),
02562 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop(),
02563 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop(),
02564 phycon.te, iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().damp(), eps );
02565 }
02566
02567 else if( strncmp( save.chSave[ipPun] , "HYD", 3 ) == 0 )
02568 {
02569
02570 if( strcmp(save.chSave[ipPun],"HYDc") == 0 )
02571 {
02572 if( strcmp(chTime,"LAST") != 0 )
02573 {
02574
02575 fprintf( save.ipPnunit[ipPun],
02576 " %.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
02577 radius.depth_mid_zone, phycon.te, dense.gas_phase[ipHYDROGEN], dense.eden,
02578 dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN],
02579 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN],
02580 hmi.H2_total/dense.gas_phase[ipHYDROGEN],
02581 findspecieslocal("H2+")->den/dense.gas_phase[ipHYDROGEN],
02582 findspecieslocal("H3+")->den/dense.gas_phase[ipHYDROGEN],
02583 findspecieslocal("H-")->den/dense.gas_phase[ipHYDROGEN] );
02584 }
02585 }
02586
02587 else if( strcmp(save.chSave[ipPun],"HYDi") == 0 )
02588 {
02589 if( strcmp(chTime,"LAST") != 0 )
02590 {
02591
02592
02593 RateInter = 0.;
02594 stage = iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ColIoniz*dense.eden*iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop();
02595 fref = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc*iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop();
02596 fout = iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop();
02597
02598 for( ion=ipH2s; ion < iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_local; ion++ )
02599 {
02600
02601 RateInter +=
02602 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Aul()*
02603 (iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Pesc() +
02604 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Pelec_esc() +
02605 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Pdest());
02606
02607 fref += iso_sp[ipH_LIKE][ipHYDROGEN].fb[ion].gamnc*iso_sp[ipH_LIKE][ipHYDROGEN].st[ion].Pop();
02608
02609 stage += iso_sp[ipH_LIKE][ipHYDROGEN].fb[ion].ColIoniz*dense.eden*
02610 iso_sp[ipH_LIKE][ipHYDROGEN].st[ion].Pop();
02611 fout += iso_sp[ipH_LIKE][ipHYDROGEN].st[ion].Pop();
02612 }
02613
02614
02615 stage /= dense.xIonDense[ipHYDROGEN][1];
02616 fref /= dense.xIonDense[ipHYDROGEN][1];
02617 fout /= dense.xIonDense[ipHYDROGEN][1];
02618
02619 fprintf( save.ipPnunit[ipPun], "hion\t%4ld\t%.2e\t%.2e\t%.2e",
02620 nzone,
02621
02622 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc,
02623 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ColIoniz* dense.EdenHCorr,
02624 ionbal.RateRecomTot[ipHYDROGEN][0] );
02625
02626 fprintf( save.ipPnunit[ipPun], "\t%.2e",
02627 iso_sp[ipH_LIKE][ipHYDROGEN].RadRec_caseB );
02628
02629 fprintf( save.ipPnunit[ipPun],
02630 "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
02631 dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0],
02632 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc/(ionbal.RateRecomTot[ipHYDROGEN][0]),
02633 iso_sp[ipH_LIKE][ipHYDROGEN].fb[1].RadRecomb[ipRecEsc],
02634 RateInter,
02635 fref/MAX2(1e-37,fout),
02636 stage/MAX2(1e-37,fout),
02637
02638 safe_div( iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc*dense.xIonDense[ipHYDROGEN][0], dense.eden*dense.xIonDense[ipHYDROGEN][1] ),
02639 secondaries.csupra[ipHYDROGEN][0]);
02640
02641 GammaPrt(iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon,rfield.nflux,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipOpac,
02642 save.ipPnunit[ipPun],iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc*
02643 0.05);
02644 }
02645 }
02646
02647 else if( strcmp(save.chSave[ipPun],"HYDp") == 0 )
02648 {
02649 if( strcmp(chTime,"LAST") != 0 )
02650 {
02651
02652
02653 fprintf( save.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e",
02654 radius.depth_mid_zone,
02655 dense.xIonDense[ipHYDROGEN][0],
02656 dense.xIonDense[ipHYDROGEN][1] );
02657
02658
02659 for( j=ipH1s; j < iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_local-1; j++ )
02660 {
02661 fprintf( save.ipPnunit[ipPun], "\t%.2e",
02662 iso_sp[ipH_LIKE][ipHYDROGEN].st[j].Pop() );
02663 }
02664 fprintf( save.ipPnunit[ipPun], "\n" );
02665 }
02666 }
02667
02668 else if( strcmp(save.chSave[ipPun],"HYDl") == 0 )
02669 {
02670 if( strcmp(chTime,"LAST") == 0 )
02671 {
02672
02673
02674 for( ipHi=1; ipHi<iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_local -
02675 iso_sp[ipH_LIKE][ipHYDROGEN].nCollapsed_local; ++ipHi )
02676 {
02677 for( ipLo=0; ipLo<ipHi; ++ipLo )
02678 {
02679 if( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipHi,ipLo).ipCont() < 0 )
02680 continue;
02681 fprintf(save.ipPnunit[ipPun], "%li\t%li\t%li\t%li\t%.4e\t%.2e\n",
02682 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipHi].n(),
02683 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipHi].l(),
02684 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipLo].n(),
02685 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipLo].l(),
02686 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipHi,ipLo).EnergyRyd(),
02687 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipHi,ipLo).Emis().TauIn() );
02688 }
02689 }
02690 }
02691 }
02692
02693
02694 else if( strcmp(save.chSave[ipPun],"HYDL") == 0 )
02695 {
02696 if( strcmp(chTime,"LAST") != 0 )
02697 {
02698
02699 double popul = iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop()/SDIV(iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop());
02700
02701 texc = TexcLine( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) );
02702 fprintf( save.ipPnunit[ipPun],
02703 "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
02704 radius.depth_mid_zone,
02705 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauIn(),
02706 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauTot(),
02707 popul,
02708 texc,
02709 phycon.te,
02710 texc/phycon.te ,
02711 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pesc(),
02712 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pdest(),
02713 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().pump(),
02714 opac.opacity_abs[iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).ipCont()-1],
02715 opac.albedo[iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).ipCont()-1] );
02716 }
02717 }
02718
02719 else if( strcmp(save.chSave[ipPun],"HYDr") == 0 )
02720 {
02721
02722 TempChange(2500.f, false);
02723 while( phycon.te <= 20000. )
02724 {
02725 double r1;
02726 double ThinCoolingCaseB;
02727
02728 r1 = HydroRecCool(1,0);
02729 ThinCoolingCaseB = pow(10.,((-25.859117 +
02730 0.16229407*phycon.telogn[0] +
02731 0.34912863*phycon.telogn[1] -
02732 0.10615964*phycon.telogn[2])/(1. +
02733 0.050866793*phycon.telogn[0] -
02734 0.014118924*phycon.telogn[1] +
02735 0.0044980897*phycon.telogn[2] +
02736 6.0969594e-5*phycon.telogn[3])))/phycon.te;
02737
02738 fprintf( save.ipPnunit[ipPun], " %10.2e\t",
02739 phycon.te);
02740 fprintf( save.ipPnunit[ipPun], " %10.2e\t",
02741 (r1+ThinCoolingCaseB)/(BOLTZMANN*phycon.te) );
02742
02743 fprintf( save.ipPnunit[ipPun], " %10.2e\t",
02744 r1/(BOLTZMANN*phycon.te));
02745
02746 fprintf( save.ipPnunit[ipPun], " %10.2e\n",
02747 ThinCoolingCaseB/(BOLTZMANN*phycon.te));
02748
02749 TempChange(phycon.te *2.f , false);
02750 }
02751
02752 fprintf(ioQQQ , "save agn now exits since solution is disturbed.\n");
02753 cdEXIT( EXIT_SUCCESS );
02754 }
02755 else
02756 TotalInsanity();
02757 }
02758
02759 else if( strcmp(save.chSave[ipPun],"IONI") == 0 )
02760 {
02761 if( strcmp(chTime,"LAST") == 0 )
02762 {
02763
02764 PrtMeanIon( 'i', false , save.ipPnunit[ipPun] );
02765 }
02766 }
02767
02768
02769 else if( strcmp(save.chSave[ipPun],"IONR") == 0 )
02770 {
02771 if( strcmp(chTime,"LAST") != 0 )
02772 {
02773
02774 nelem = (long)save.punarg[ipPun][0];
02775 fprintf( save.ipPnunit[ipPun],
02776 "%.5e\t%.4e\t%.4e",
02777 radius.depth_mid_zone,
02778 dense.eden ,
02779 dynamics.Rate);
02780
02781
02782 for( ion=0; ion<nelem+1; ++ion )
02783 {
02784 fprintf( save.ipPnunit[ipPun],
02785 "\t%.4e\t%.4e\t%.4e\t%.4e",
02786 dense.xIonDense[nelem][ion] ,
02787 ionbal.RateIonizTot(nelem,ion) ,
02788 ionbal.RateRecomTot[nelem][ion] ,
02789 dynamics.Source[nelem][ion] );
02790 }
02791 fprintf( save.ipPnunit[ipPun], "\n");
02792 }
02793 }
02794
02795 else if( strcmp(save.chSave[ipPun]," IP ") == 0 )
02796 {
02797 if( strcmp(chTime,"LAST") == 0 )
02798 {
02799
02800 for( nelem=0; nelem < LIMELM; nelem++ )
02801 {
02802 int ion_big;
02803 double energy;
02804
02805
02806 const int NELEM_LINE = 10;
02807
02808 for( ion_big=0; ion_big<=nelem; ion_big += NELEM_LINE )
02809 {
02810 int ion_limit = MIN2(ion_big+NELEM_LINE-1,nelem);
02811
02812
02813 fprintf( save.ipPnunit[ipPun],
02814 "\n%2.2s", elementnames.chElementSym[nelem]);
02815
02816
02817 for( ion=ion_big; ion <= ion_limit; ++ion )
02818 {
02819 fprintf( save.ipPnunit[ipPun], "\t%4ld", ion+1 );
02820 }
02821 fprintf( save.ipPnunit[ipPun], "\n" );
02822
02823
02824 ASSERT( ion_limit < LIMELM );
02825
02826 for( ips=0; ips < Heavy.nsShells[nelem][ion_big]; ips++ )
02827 {
02828
02829
02830 fprintf( save.ipPnunit[ipPun], "%2.2s", Heavy.chShell[ips]);
02831
02832
02833 for( ion=ion_big; ion<=ion_limit; ++ion )
02834 {
02835
02836
02837
02838 if( ips >= Heavy.nsShells[nelem][ion] )
02839 break;
02840
02841
02842 energy = t_ADfA::Inst().ph1(ips,nelem-ion,nelem,0);
02843
02844
02845 if( energy < 10. )
02846 {
02847 fprintf( save.ipPnunit[ipPun], "\t%6.3f", energy );
02848 }
02849 else if( energy < 100. )
02850 {
02851 fprintf( save.ipPnunit[ipPun], "\t%6.2f", energy );
02852 }
02853 else if( energy < 1000. )
02854 {
02855 fprintf( save.ipPnunit[ipPun], "\t%6.1f", energy );
02856 }
02857 else
02858 {
02859 fprintf( save.ipPnunit[ipPun], "\t%6ld", (long)(energy) );
02860 }
02861 }
02862
02863
02864 fprintf( save.ipPnunit[ipPun], "\n" );
02865 }
02866 }
02867 }
02868 }
02869 }
02870
02871 else if( strcmp(save.chSave[ipPun],"LINC") == 0 )
02872 {
02873
02874 if( strcmp(chTime,"LAST") != 0 )
02875 {
02876 save_line(save.ipPnunit[ipPun],"PUNC",
02877 save.lgEmergent[ipPun]);
02878 }
02879 }
02880
02881 else if( strcmp(save.chSave[ipPun],"LIND") == 0 )
02882 {
02883
02884 SaveLineData(save.ipPnunit[ipPun]);
02885 }
02886
02887 else if( strcmp(save.chSave[ipPun],"LINL") == 0 )
02888 {
02889
02890 bool lgPrintAll=false;
02891
02892 if( save.punarg[ipPun][0]>0. )
02893 lgPrintAll = true;
02894 prt_LineLabels(save.ipPnunit[ipPun] , lgPrintAll );
02895 }
02896
02897 else if( strcmp(save.chSave[ipPun],"LINO") == 0 )
02898 {
02899 if( strcmp(chTime,"LAST") == 0 )
02900 {
02901
02902 SaveLineStuff(save.ipPnunit[ipPun],"optical" , save.punarg[ipPun][0]);
02903 }
02904 }
02905
02906 else if( strcmp(save.chSave[ipPun],"LINP") == 0 )
02907 {
02908 if( strcmp(chTime,"LAST") != 0 )
02909 {
02910 static bool lgFirst=true;
02911
02912
02913
02914 SaveLineStuff(save.ipPnunit[ipPun],"populat" , save.punarg[ipPun][0]);
02915 if( lgFirst )
02916 {
02917 lgFirst = false;
02918 SaveLineStuff(save.ipPnunit[ipPun],"populat" , save.punarg[ipPun][0]);
02919 }
02920 }
02921 }
02922
02923 else if( strcmp(save.chSave[ipPun],"LINS") == 0 )
02924 {
02925
02926 if( strcmp(chTime,"LAST") != 0 )
02927 {
02928 save_line(save.ipPnunit[ipPun],"PUNS",
02929 save.lgEmergent[ipPun]);
02930 }
02931 }
02932
02933 else if( strcmp(save.chSave[ipPun],"LINR") == 0 )
02934 {
02935
02936 if( strcmp(chTime,"LAST") != 0 )
02937 Save_Line_RT( save.ipPnunit[ipPun]);
02938 }
02939
02940 else if( strcmp(save.chSave[ipPun],"LINA") == 0 )
02941 {
02942
02943 if( strcmp(chTime,"LAST") == 0 )
02944 {
02945
02946 for( j=0; j < LineSave.nsum; j++ )
02947 {
02948 if( LineSv[j].wavelength > 0. &&
02949 LineSv[j].SumLine[0] > 0. )
02950 {
02951
02952 fprintf( save.ipPnunit[ipPun], "%12.5e",
02953 AnuUnit((realnum)RYDLAM/LineSv[j].wavelength) );
02954
02955 fprintf( save.ipPnunit[ipPun], "\t%4.4s ",
02956 LineSv[j].chALab );
02957
02958 prt_wl( save.ipPnunit[ipPun], LineSv[j].wavelength );
02959
02960 fprintf( save.ipPnunit[ipPun], "\t%8.3f",
02961 log10(SDIV(LineSv[j].SumLine[0]) ) + radius.Conv2PrtInten );
02962
02963 fprintf( save.ipPnunit[ipPun], "\t%8.3f",
02964 log10(SDIV(LineSv[j].SumLine[1]) ) + radius.Conv2PrtInten );
02965
02966 fprintf( save.ipPnunit[ipPun], " \t%c\n",
02967 LineSv[j].chSumTyp);
02968 }
02969 }
02970 }
02971 }
02972
02973 else if( strcmp(save.chSave[ipPun],"LINI") == 0 )
02974 {
02975 if( strcmp(chTime,"LAST") == 0 &&
02976 (nzone/save.LinEvery)*save.LinEvery != nzone )
02977 {
02978
02979
02980 SaveLineIntensity(save.ipPnunit[ipPun] , ipPun , save.punarg[ipPun][0] );
02981 }
02982 else if( strcmp(chTime,"LAST") != 0 )
02983 {
02984
02985 if( (save.lgLinEvery && nzone == 1) ||
02986 (nzone/save.LinEvery)*save.LinEvery == nzone )
02987 {
02988
02989
02990 SaveLineIntensity(save.ipPnunit[ipPun] , ipPun , save.punarg[ipPun][0]);
02991 }
02992 }
02993 }
02994
02995 else if( strcmp( save.chSave[ipPun],"LEIL") == 0)
02996 {
02997
02998
02999 if( strcmp(chTime,"LAST") == 0 )
03000 {
03001 double absval , rel;
03002 long int n;
03003
03004
03005
03006
03007 const int NLINE_H2 = 31;
03008
03009 const int NLINE_NOTH_H2 = 5;
03010
03011 char chLabel[NLINE_NOTH_H2][5]=
03012 {"C 2", "O 1", "O 1","C 1", "C 1" };
03013 double Wl[NLINE_NOTH_H2]=
03014 {157.6 , 63.17 , 145.5 ,609.2 , 369.7 , };
03015
03016
03017
03018 double Wl_H2[NLINE_H2]=
03019 {2.121 ,
03020 28.213, 17.03 , 12.28 , 9.662 , 8.024 , 6.907 , 6.107 , 5.510 , 5.051 , 4.693 ,
03021 4.408 , 4.180 , 3.996 , 3.845 , 3.723 , 3.625 , 3.546 , 3.485 , 3.437 , 3.403 ,
03022 3.380 , 3.368 , 3.365 , 3.371 , 3.387 , 3.410 , 3.441 , 3.485 , 3.542 , 3.603};
03023
03024 for( n=0; n<NLINE_NOTH_H2; ++n )
03025 {
03026 fprintf(save.ipPnunit[ipPun], "%s\t%.2f",chLabel[n] , Wl[n]);
03027
03028
03029 if( cdLine( chLabel[n] , (realnum)(Wl[n]*1e4) , &absval , &rel ) <= 0 )
03030 {
03031 fprintf(save.ipPnunit[ipPun], " did not find\n");
03032 }
03033 else
03034 {
03035 fprintf(save.ipPnunit[ipPun], "\t%.3e\t%.3e\n",pow(10.,rel),absval);
03036 }
03037 }
03038 fprintf(save.ipPnunit[ipPun], "\n\n\n");
03039
03040
03041 if( h2.lgEnabled )
03042 {
03043 fprintf(save.ipPnunit[ipPun],
03044 "Here are some of the H2 Intensities, The first one is the\n"
03045 "1-0 S(0) line and the following ones are the 0-0 S(X)\n"
03046 "lines where X goes from 0 to 29\n\n");
03047 for( n=0; n<NLINE_H2; ++n )
03048 {
03049 fprintf(save.ipPnunit[ipPun], "%s\t%.3f","H2 " , Wl_H2[n]);
03050
03051 if( cdLine( "H2 " , (realnum)(Wl_H2[n]*1e4) , &absval , &rel ) <= 0 )
03052 {
03053 fprintf(save.ipPnunit[ipPun], " did not find\n");
03054 }
03055 else
03056 {
03057 fprintf(save.ipPnunit[ipPun], "\t%.3e\t%.3e\n",pow(10.,rel),absval);
03058 }
03059 }
03060 }
03061 }
03062 }
03063
03064 else if( strcmp( save.chSave[ipPun],"LEIS") == 0)
03065 {
03066 if( strcmp(chTime,"LAST") != 0 )
03067 {
03068
03069 double col_ci , col_oi , col_cii, col_heii;
03070 if( cdColm("carb" , 1 , &col_ci ) )
03071 TotalInsanity();
03072 if( cdColm("carb" , 2 , &col_cii ) )
03073 TotalInsanity();
03074 if( cdColm("oxyg" , 1 , &col_oi ) )
03075 TotalInsanity();
03076 if( cdColm("heli" , 2 , &col_heii ) )
03077 TotalInsanity();
03078
03079 fprintf( save.ipPnunit[ipPun],
03080 "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
03081 "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
03082 "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
03083 "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
03084 "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
03085
03086 radius.depth_mid_zone,
03087
03088 0.00,
03089
03090 rfield.extin_mag_V_point,
03091
03092 phycon.te ,
03093 dense.xIonDense[ipHYDROGEN][0],
03094 hmi.H2_total,
03095 dense.xIonDense[ipCARBON][0],
03096 dense.xIonDense[ipCARBON][1],
03097 dense.xIonDense[ipOXYGEN][0],
03098 findspecieslocal("CO")->den,
03099 findspecieslocal("O2")->den,
03100 findspecieslocal("CH")->den,
03101 findspecieslocal("OH")->den,
03102 dense.eden,
03103 dense.xIonDense[ipHELIUM][1],
03104 dense.xIonDense[ipHYDROGEN][1],
03105 findspecieslocal("H3+")->den,
03106 colden.colden[ipCOL_H0],
03107 colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s],
03108 col_ci,
03109 col_cii,
03110 col_oi,
03111 findspecieslocal("CO")->column,
03112 findspecieslocal("O2")->column,
03113 findspecieslocal("CH")->column,
03114 findspecieslocal("OH")->column,
03115 colden.colden[ipCOL_elec],
03116 col_heii,
03117 colden.colden[ipCOL_Hp],
03118 colden.colden[ipCOL_H3p],
03119 hmi.H2_Solomon_dissoc_rate_used_H2g ,
03120 gv.rate_h2_form_grains_used_total,
03121 hmi.H2_photodissoc_used_H2g,
03122 hmi.UV_Cont_rel2_Draine_DB96_depth,
03123
03124 mole.findrk("PHOTON,CO=>C,O"),
03125
03126 ionbal.PhotoRate_Shell[ipCARBON][0][2][0],
03127
03128 thermal.htot,
03129
03130 thermal.ctot,
03131
03132 thermal.heating[0][13],
03133
03134 MAX2(0.,gv.GasCoolColl),
03135
03136 -1.*MIN2(0.,gv.GasCoolColl),
03137
03138 thermal.heating[0][9],
03139
03140 hmi.HeatH2Dish_used,
03141
03142 hmi.HeatH2Dexc_used ,
03143
03144 thermal.heating[0][24] ,
03145
03146 thermal.heating[1][6] ,
03147
03148 thermal.heating[ipMAGNESIUM][0],
03149 thermal.heating[ipSULPHUR][0],
03150 thermal.heating[ipSILICON][0],
03151 thermal.heating[ipIRON][0],
03152 thermal.heating[ipSODIUM][0],
03153 thermal.heating[ipALUMINIUM][0],
03154 thermal.heating[ipCARBON][0],
03155 TauLines[ipT610].Coll().cool(),
03156 TauLines[ipT370].Coll().cool(),
03157 TauLines[ipT157].Coll().cool(),
03158 TauLines[ipT63].Coll().cool(),
03159 TauLines[ipT146].Coll().cool() );
03160 }
03161 }
03162
03163 # ifdef USE_NLTE7
03164 else if( strcmp( save.chSave[ipPun],"NLTE") == 0)
03165 {
03166 if( strcmp(chTime,"LAST") == 0 )
03167 {
03168
03169 runNLTE(ipPun);
03170 }
03171 }
03172 # endif
03173
03174
03175
03176 else if( strcmp( save.chSave[ipPun],"LY1") == 0)
03177 {
03178 static bool runonce = true;
03179 if (runonce )
03180 {
03181 for( long ipHi=0; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
03182 {
03183 ASSERT(FeII.FeIINRGs[ipHi] >= 0.);
03184 ASSERT(FeII.FeIISTWT[ipHi] >= 0.);
03185 fprintf(save.ipPnunit[ipPun], "%li\t%10.3f\t%3.1f\n",ipHi+1,FeII.FeIINRGs[ipHi],FeII.FeIISTWT[ipHi]);
03186 }
03187 runonce = false;
03188 fprintf(save.ipPnunit[ipPun], "-1\n");
03189 }
03190 }
03191
03192
03193 else if( strcmp( save.chSave[ipPun],"LY2") == 0)
03194 {
03195 static bool runonce = true;
03196 if (runonce )
03197 {
03198 for( long ipLo = 0; ipLo < FeII.nFeIILevel_malloc; ipLo++)
03199 {
03200 for( long ipHi=0; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
03201 {
03202 if( FeII.FeIIAul[ipHi][ipLo] >= 0. && FeII.FeIIAul[ipHi][ipLo] != 1e-5f && FeII.FeIIAul[ipHi][ipLo] != 1e-6f )
03203 {
03204 fprintf(save.ipPnunit[ipPun], "A\t%li\t%li\t%8.2e\n",
03205 ipLo+1,ipHi+1,FeII.FeIIAul[ipHi][ipLo]);
03206 }
03207 }
03208 }
03209 runonce = false;
03210 fprintf(save.ipPnunit[ipPun], "-1\n");
03211 }
03212 }
03213
03214
03215 else if( strcmp( save.chSave[ipPun],"LY3") == 0)
03216 {
03217 static realnum tt[8]={1e3f,3e3f,5e3f,7e3f,1e4f,12e3f,15e3f,2e4f};
03218 static bool runonce = true;
03219 if (runonce )
03220 {
03221 fprintf(save.ipPnunit[ipPun], "TEMP\t8");
03222 for( int k = 0; k < 8; k++)
03223 {
03224 fprintf(save.ipPnunit[ipPun], "\t%8.2e",tt[k]);
03225 }
03226 fprintf(save.ipPnunit[ipPun], "\n");
03227 for( long ipHi = 1; ipHi < FeII.nFeIILevel_malloc; ipHi++)
03228 {
03229 for( long ipLo=0; ipLo < ipHi; ipLo++ )
03230 {
03231 bool skipLine = false;
03232 for( int k = 0; k < 8; k++)
03233 {
03234 if( FeII.FeIIColl[ipHi][ipLo][k] == -2. || FeII.FeIIColl[ipHi][ipLo][k] == 0.)
03235 {
03236 skipLine = true;
03237 break;
03238 }
03239 }
03240 if( !skipLine )
03241 {
03242 fprintf(save.ipPnunit[ipPun], "CSELECTRON\t%li\t%li",ipLo+1,ipHi+1);
03243 for( int k = 0; k < 8; k++)
03244 {
03245 fprintf(save.ipPnunit[ipPun], "\t%8.2e",FeII.FeIIColl[ipHi][ipLo][k]);
03246 }
03247 fprintf(save.ipPnunit[ipPun], "\n");
03248 }
03249 }
03250 }
03251 runonce = false;
03252 fprintf(save.ipPnunit[ipPun], "-1\n");
03253 }
03254 }
03255
03256 else if( strcmp( save.chSave[ipPun],"LLST") == 0)
03257 {
03258
03259 if( strcmp(chTime,"LAST") == 0 )
03260 {
03261 fprintf( save.ipPnunit[ipPun], "iteration %li" , iteration );
03262 if( save.punarg[ipPun][1] )
03263 fprintf( save.ipPnunit[ipPun], "\n" );
03264
03265
03266 if( save.nLineList[ipPun] < 0 )
03267 TotalInsanity();
03268
03269 int LineType = 0;
03270 if( save.lgEmergent[ipPun] )
03271 LineType = 1;
03272 if( save.lgCumulative[ipPun] )
03273 LineType += 2;
03274
03275
03276 for( j=0; j<save.nLineList[ipPun]; ++j )
03277 {
03278 double relative , absolute, PrtQuantity;
03279 if( (cdLine( save.chLineListLabel[ipPun][j] ,
03280 save.wlLineList[ipPun][j] ,
03281 &relative , &absolute , LineType ) ) <=0 )
03282 {
03283 if( !h2.lgEnabled && strncmp( save.chLineListLabel[ipPun][j] , "H2 " , 4 )==0 )
03284 {
03285 static bool lgMustPrintFirstTime = true;
03286 if( lgMustPrintFirstTime )
03287 {
03288
03289 fprintf( ioQQQ,"Did not find an H2 line, the large model is not "
03290 "included, so I will ignore it. Log intensity set to -30.\n" );
03291 fprintf( ioQQQ,"I will totally ignore any future missed H2 lines\n");
03292 lgMustPrintFirstTime = false;
03293 }
03294 relative = -30.f;
03295 absolute = -30.f;
03296 }
03297 else if( lgAbort )
03298 {
03299
03300 relative = -30.f;
03301 absolute = -30.f;
03302 }
03303 else
03304 {
03305 fprintf(ioQQQ,"DISASTER - did not find a line in the Line List table\n");
03306 cdEXIT(EXIT_FAILURE);
03307 }
03308 }
03309
03310
03311
03312
03313
03314 if( save.punarg[ipPun][0] > 0 )
03315 PrtQuantity = pow(10. , absolute);
03316 else
03317 PrtQuantity = relative;
03318
03319
03320 if( save.punarg[ipPun][1] )
03321 {
03322
03323 if( save.lgLineListRatio[ipPun] && is_odd(j) )
03324 fprintf( save.ipPnunit[ipPun] , "/" );
03325
03326 fprintf( save.ipPnunit[ipPun], "%s ", save.chLineListLabel[ipPun][j] );
03327 char chTemp[MAX_HEADER_SIZE];
03328 sprt_wl( chTemp, save.wlLineList[ipPun][j] );
03329 fprintf( save.ipPnunit[ipPun], "%s ", chTemp );
03330 }
03331
03332
03333
03334 if( save.lgLineListRatio[ipPun] )
03335 {
03336
03337 static double SaveQuantity = 0;
03338 if( is_odd(j) )
03339 fprintf( save.ipPnunit[ipPun], "\t%.4e" ,
03340 SaveQuantity / SDIV( PrtQuantity ) );
03341 else
03342 SaveQuantity = PrtQuantity;
03343 }
03344 else
03345 {
03346 fprintf( save.ipPnunit[ipPun], "\t%.4e" , PrtQuantity );
03347 }
03348
03349 if( save.punarg[ipPun][1] )
03350 {
03351 if( !save.lgLineListRatio[ipPun] ||
03352 is_odd(j) )
03353 fprintf( save.ipPnunit[ipPun], "\n" );
03354 }
03355 }
03356 fprintf( save.ipPnunit[ipPun], "\n" );
03357 }
03358 }
03359
03360 else if( strcmp( save.chSave[ipPun],"CHRT") == 0)
03361 {
03362
03363 if( strcmp(chTime,"LAST") != 0 )
03364 {
03365 bool lgHeader, lgData;
03366 if( save.lgPunHeader[ipPun] )
03367 {
03368 lgHeader = true;
03369 lgData = false;
03370 mole_punch(save.ipPnunit[ipPun],save.optname[ipPun].c_str(),save.chSaveArgs[ipPun],lgHeader,lgData,radius.depth_mid_zone);
03371 }
03372 save.lgPunHeader[ipPun] = false;
03373 lgHeader = false;
03374 lgData = true;
03375 mole_punch(save.ipPnunit[ipPun],save.optname[ipPun].c_str(),save.chSaveArgs[ipPun],lgHeader,lgData,radius.depth_mid_zone);
03376 }
03377 }
03378
03379 else if( strcmp(save.chSave[ipPun],"MAP ") == 0 )
03380 {
03381
03382
03383 if( !hcmap.lgMapDone &&
03384 (nzone == hcmap.MapZone || strcmp(chTime,"LAST") == 0 ) )
03385 {
03386 lgTlkSav = called.lgTalk;
03387 called.lgTalk = cpu.i().lgMPI_talk();
03388 hcmap.lgMapBeingDone = true;
03389 map_do(save.ipPnunit[ipPun] , " map");
03390 called.lgTalk = lgTlkSav;
03391 }
03392 }
03393
03394 else if( strcmp(save.chSave[ipPun],"MOLE") == 0 )
03395 {
03396 if( save.lgPunHeader[ipPun] )
03397 {
03398 fprintf( save.ipPnunit[ipPun],
03399 "#molecular species will follow:\n");
03400 fprintf( save.ipPnunit[ipPun],
03401 "#depth\tAV(point)\tAV(extend)\tCO diss rate\tC recom rate");
03402
03403 for(i=0; i<mole_global.num_calc; ++i )
03404 {
03405 fprintf( save.ipPnunit[ipPun], "\t%s", mole_global.list[i]->label.c_str() );
03406 }
03407 fprintf ( save.ipPnunit[ipPun], "\n");
03408 save.lgPunHeader[ipPun] = false;
03409 }
03410 if( strcmp(chTime,"LAST") != 0 )
03411 {
03412
03413 fprintf( save.ipPnunit[ipPun], "%.5e\t" , radius.depth_mid_zone );
03414
03415
03416 fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_point);
03417
03418
03419 fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended);
03420
03421
03422 fprintf( save.ipPnunit[ipPun], "%.5e\t" , mole.findrk("PHOTON,CO=>C,O") );
03423
03424
03425 fprintf( save.ipPnunit[ipPun], "%.5e" , ionbal.RateRecomTot[ipCARBON][0] );
03426
03427
03428 for(j=0; j<mole_global.num_calc; ++j )
03429 {
03430 fprintf(save.ipPnunit[ipPun],"\t%.2e",mole.species[j].den );
03431 }
03432
03433 fprintf(save.ipPnunit[ipPun],"\n");
03434 }
03435 }
03436
03437 else if( strcmp(save.chSave[ipPun],"OPAC") == 0 )
03438 {
03439
03440 if( save.lgSaveEveryZone[ipPun] || strcmp(chTime,"LAST") == 0 )
03441 save_opacity(save.ipPnunit[ipPun],ipPun);
03442 }
03443
03444
03445 else if( strcmp(save.chSave[ipPun],"OPTc") == 0 )
03446 {
03447 if( save.lgSaveEveryZone[ipPun] || strcmp(chTime,"LAST") == 0 )
03448 {
03449 for( j=0; j < rfield.nflux; j++ )
03450 {
03451 fprintf( save.ipPnunit[ipPun],
03452 "%13.5e\t%.3e\t%12.4e\t%.3e\n",
03453 AnuUnit(rfield.AnuOrg[j]),
03454 opac.TauAbsFace[j]+opac.TauScatFace[j],
03455 opac.TauAbsFace[j],
03456 opac.TauScatFace[j] );
03457 }
03458 }
03459 }
03460
03461
03462 else if( strcmp(save.chSave[ipPun],"OPTf") == 0 )
03463 {
03464 if( save.lgSaveEveryZone[ipPun] || strcmp(chTime,"LAST") == 0 )
03465 {
03466 long nu_hi , nskip;
03467 if( save.punarg[ipPun][0] > 0. )
03468
03469 j = ipFineCont( save.punarg[ipPun][0] );
03470 else
03471 j = 0;
03472
03473
03474 if( save.punarg[ipPun][1]> 0. )
03475 nu_hi = ipFineCont( save.punarg[ipPun][1]);
03476 else
03477 nu_hi = rfield.nfine;
03478
03479
03480
03481 nskip = (long)save.punarg[ipPun][2];
03482 do
03483 {
03484 realnum sum1 = rfield.fine_opt_depth[j];
03485 realnum sum2 = rfield.fine_opac_zone[j];
03486
03487 realnum xnu = rfield.fine_anu[j];
03488 for( jj=1; jj<nskip; ++jj )
03489 {
03490 sum1 += rfield.fine_opt_depth[j+jj];
03491 sum2 += rfield.fine_opac_zone[j+jj];
03492 xnu += rfield.fine_anu[j+jj];
03493 }
03494 if( sum2 > 0. )
03495 fprintf( save.ipPnunit[ipPun],
03496 "%12.6e\t%.3e\t%.3e\n",
03497 AnuUnit(xnu/nskip),
03498 sum1/nskip ,
03499 sum2/nskip);
03500 j += nskip;
03501 }while( j < nu_hi );
03502 }
03503 }
03504
03505 else if( strcmp(save.chSave[ipPun]," OTS") == 0 )
03506 {
03507 ConMax = 0.;
03508 xLinMax = 0.;
03509 opConSum = 0.;
03510 opLinSum = 0.;
03511 ipLinMax = 1;
03512 ipConMax = 1;
03513
03514 for( j=0; j < rfield.nflux; j++ )
03515 {
03516 opConSum += rfield.otscon[j]*opac.opacity_abs[j];
03517 opLinSum += rfield.otslin[j]*opac.opacity_abs[j];
03518 if( rfield.otslin[j]*opac.opacity_abs[j] > xLinMax )
03519 {
03520 xLinMax = rfield.otslin[j]*opac.opacity_abs[j];
03521 ipLinMax = j+1;
03522 }
03523 if( rfield.otscon[j]*opac.opacity_abs[j] > ConMax )
03524 {
03525 ConMax = rfield.otscon[j]*opac.opacity_abs[j];
03526 ipConMax = j+1;
03527 }
03528 }
03529 fprintf( save.ipPnunit[ipPun],
03530 "tot con lin=%.2e%.2e lin=%.4s%.4e%.2e con=%.4s%.4e%.2e\n",
03531 opConSum, opLinSum, rfield.chLineLabel[ipLinMax-1]
03532 , rfield.anu[ipLinMax-1], xLinMax, rfield.chContLabel[ipConMax-1]
03533 , rfield.anu[ipConMax-1], ConMax );
03534 }
03535
03536 else if( strcmp(save.chSave[ipPun],"OVER") == 0 )
03537 {
03538
03539
03540 double toosmall = SMALLFLOAT ,
03541 hold;
03542
03543
03544
03545 if( strcmp(chTime,"LAST") != 0 )
03546 {
03547
03548
03549 fprintf( save.ipPnunit[ipPun], "%.5e\t", radius.depth_mid_zone );
03550
03551
03552 if(dynamics.Cool() > dynamics.Heat())
03553 {
03554 fprintf( save.ipPnunit[ipPun], "%.4f\t%.3f",
03555 log10(phycon.te), log10(SDIV(thermal.htot-dynamics.Heat())) );
03556 }
03557 else
03558 {
03559 double diff = fabs(thermal.htot-dynamics.Cool());
03560 fprintf( save.ipPnunit[ipPun], "%.4f\t%.3f",
03561 log10(phycon.te), log10(SDIV(diff)) );
03562 }
03563
03564
03565 fprintf( save.ipPnunit[ipPun], "\t%.4f\t%.4f",
03566 log10(dense.gas_phase[ipHYDROGEN]), log10(dense.eden) );
03567
03568
03569 fprintf( save.ipPnunit[ipPun], "\t%.4f",
03570
03571 log10(MAX2(toosmall,2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN])) );
03572
03573
03574 fprintf( save.ipPnunit[ipPun], "\t%.4f\t%.4f",
03575 log10(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN])),
03576 log10(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN])) );
03577
03578
03579 for( j=1; j <= 3; j++ )
03580 {
03581 double arg1 = SDIV(dense.gas_phase[ipHELIUM]);
03582 arg1 = MAX2(toosmall,dense.xIonDense[ipHELIUM][j-1]/arg1 );
03583 fprintf( save.ipPnunit[ipPun], "\t%.4f",
03584 log10(arg1) );
03585 }
03586
03587
03588 hold = SDIV(dense.gas_phase[ipCARBON]);
03589 hold = findspecieslocal("CO")->den/hold;
03590 hold = MAX2(toosmall, hold );
03591 fprintf( save.ipPnunit[ipPun], "\t%.4f", log10(hold) );
03592
03593
03594 for( j=1; j <= 4; j++ )
03595 {
03596 hold = SDIV(dense.gas_phase[ipCARBON]);
03597 hold = MAX2(toosmall,dense.xIonDense[ipCARBON][j-1]/hold);
03598 fprintf( save.ipPnunit[ipPun], "\t%.4f",
03599 log10(hold) );
03600 }
03601
03602
03603 for( j=1; j <= 6; j++ )
03604 {
03605 hold = SDIV(dense.gas_phase[ipOXYGEN]);
03606 hold = MAX2(toosmall,dense.xIonDense[ipOXYGEN][j-1]/hold);
03607 fprintf( save.ipPnunit[ipPun], "\t%.4f",
03608 log10(hold) );
03609 }
03610
03611
03612 hold = SDIV(dense.gas_phase[ipOXYGEN]);
03613 hold = findspecieslocal("H2O")->den/hold;
03614 hold = MAX2(toosmall, hold );
03615 fprintf( save.ipPnunit[ipPun], "\t%.4f", log10(hold) );
03616
03617
03618 fprintf( save.ipPnunit[ipPun], "\t%.2e" , rfield.extin_mag_V_point);
03619
03620
03621 fprintf( save.ipPnunit[ipPun], "\t%.2e\n" , rfield.extin_mag_V_extended);
03622 }
03623 }
03624
03625 else if( strcmp(save.chSave[ipPun]," PDR") == 0 )
03626 {
03627
03628 if( strcmp(chTime,"LAST") != 0 )
03629 {
03630
03631
03632
03633
03634
03635
03636
03637 fprintf( save.ipPnunit[ipPun],
03638 "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t",
03639 radius.depth_mid_zone,
03640
03641 colden.colden[ipCOL_HTOT],
03642 phycon.te,
03643
03644 dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN],
03645
03646 2.*findspecieslocal("H2")->den/dense.gas_phase[ipHYDROGEN],
03647 2.*findspecieslocal("H2*")->den/dense.gas_phase[ipHYDROGEN],
03648
03649 dense.xIonDense[ipCARBON][0]/SDIV(dense.gas_phase[ipCARBON]),
03650 findspecieslocal("CO")->den/SDIV(dense.gas_phase[ipCARBON]),
03651 findspecieslocal("H2O")->den/SDIV(dense.gas_phase[ipOXYGEN]),
03652
03653 hmi.UV_Cont_rel2_Habing_TH85_depth);
03654
03655
03656 fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_point);
03657
03658
03659 fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended);
03660
03661
03662 fprintf( save.ipPnunit[ipPun], "%.2e\n" , opac.TauAbsGeo[0][rfield.ipV_filter] );
03663 }
03664 }
03665
03666
03667 else if( strcmp(save.chSave[ipPun],"PERF") == 0 )
03668 {
03669 if( strcmp(chTime,"LAST") != 0 )
03670 {
03671 static double ElapsedTime , ZoneTime;
03672 if( nzone<=1 )
03673 {
03674 ElapsedTime = cdExecTime();
03675 ZoneTime = 0.;
03676 }
03677 else
03678 {
03679 double t = cdExecTime();
03680 ZoneTime = t - ElapsedTime;
03681 ElapsedTime = t;
03682 }
03683
03684
03685 fprintf( save.ipPnunit[ipPun], " %ld\t%.3f\t%.2f\t%li",
03686 nzone, ZoneTime , ElapsedTime, conv.nPres2Ioniz );
03687
03688 for( long i=0; i<NTYPES; ++i )
03689 fprintf( save.ipPnunit[ipPun], "\t%li", conv.getCounterZone(i) );
03690 fprintf( save.ipPnunit[ipPun], "\n" );
03691 }
03692 }
03693
03694 else if( strcmp(save.chSave[ipPun],"PHYS") == 0 )
03695 {
03696 if( strcmp(chTime,"LAST") != 0 )
03697 {
03698
03699 fprintf( save.ipPnunit[ipPun], "%.5e\t%.4e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
03700 radius.depth_mid_zone, phycon.te, dense.gas_phase[ipHYDROGEN],
03701 dense.eden, thermal.htot, wind.AccelTotalOutward, geometry.FillFac );
03702 }
03703 }
03704
03705 else if( strcmp(save.chSave[ipPun],"PRES") == 0 )
03706 {
03707
03708 if( strcmp(chTime,"LAST") != 0 )
03709 {
03710 fprintf( save.ipPnunit[ipPun],
03711 "%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%c\n",
03712
03713 radius.depth_mid_zone,
03714
03715 pressure.PresTotlError*100.,
03716
03717 pressure.PresTotlCurr,
03718
03719
03720
03721 pressure.PresTotlInit + pressure.PresInteg - pressure.pinzon,
03722
03723 pressure.PresTotlInit,
03724
03725 pressure.PresGasCurr,
03726
03727 pressure.PresRamCurr,
03728
03729 pressure.pres_radiation_lines_curr,
03730
03731 pressure.PresInteg - pressure.pinzon,
03732
03733 wind.windv/1e5,
03734
03735 timesc.sound_speed_adiabatic/1e5,
03736
03737 magnetic.pressure ,
03738
03739 DoppVel.TurbVel/1e5 ,
03740
03741 pressure.PresTurbCurr*DoppVel.lgTurb_pressure ,
03742
03743 pressure.IntegRhoGravity,
03744
03745
03746
03747
03748 pressure.PresIntegElecThin-pressure.pinzon_PresIntegElecThin,
03749
03750 TorF(conv.lgConvPres) );
03751 }
03752 }
03753 else if( strcmp(save.chSave[ipPun],"PREL") == 0 )
03754 {
03755
03756 fprintf( save.ipPnunit[ipPun],
03757 "%.5e\t%.3e\t%.3e\t",
03758
03759 radius.depth_mid_zone ,
03760 pressure.PresTotlCurr,
03761 pressure.pres_radiation_lines_curr/SDIV(pressure.PresTotlCurr) );
03762 PrtLinePres(save.ipPnunit[ipPun]);
03763
03764 }
03765
03766 else if( save.chSave[ipPun][0]=='R' )
03767 {
03768
03769 if( strcmp(save.chSave[ipPun],"RADI") == 0 )
03770 {
03771
03772 if( strcmp(chTime,"LAST") != 0 )
03773 {
03774 fprintf( save.ipPnunit[ipPun], "%ld\t%.5e\t%.4e\t%.4e\n",
03775 nzone, radius.Radius_mid_zone, radius.depth_mid_zone,
03776 radius.drad );
03777 }
03778 }
03779
03780 else if( strcmp(save.chSave[ipPun],"RADO") == 0 )
03781 {
03782
03783 if( strcmp(chTime,"LAST") == 0 )
03784 {
03785 fprintf( save.ipPnunit[ipPun], "%ld\t%.5e\t%.4e\t%.4e\n",
03786 nzone, radius.Radius_mid_zone, radius.depth_mid_zone,
03787 radius.drad );
03788 }
03789 }
03790
03791 else if( strcmp(save.chSave[ipPun],"RESU") == 0 )
03792 {
03793
03794 if( strcmp(chTime,"LAST") == 0 )
03795 SaveResults(save.ipPnunit[ipPun]);
03796 }
03797 else
03798 {
03799
03800 TotalInsanity();
03801 }
03802 }
03803
03804 else if( strcmp(save.chSave[ipPun],"SECO") == 0 )
03805 {
03806
03807 if( strcmp(chTime,"LAST") != 0 )
03808 fprintf(save.ipPnunit[ipPun],
03809 "%.5e\t%.3e\t%.3e\t%.3e\n",
03810 radius.depth ,
03811 secondaries.csupra[ipHYDROGEN][0],
03812 secondaries.csupra[ipHYDROGEN][0]*2.02,
03813 secondaries.x12tot );
03814 }
03815
03816 else if( strcmp(save.chSave[ipPun],"SOUS") == 0 )
03817 {
03818
03819
03820 if( strcmp(chTime,"LAST") != 0 )
03821 {
03822 limit = MIN2(rfield.ipMaxBolt,rfield.nflux);
03823 for( j=0; j < limit; j++ )
03824 {
03825 fprintf( save.ipPnunit[ipPun],
03826 "%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
03827 AnuUnit(rfield.AnuOrg[j]),
03828 rfield.ConEmitLocal[nzone][j]/rfield.widflx[j],
03829 opac.opacity_abs[j],
03830 rfield.ConSourceFcnLocal[nzone][j],
03831 rfield.ConSourceFcnLocal[nzone][j]/plankf(j) ,
03832 safe_div(rfield.ConSourceFcnLocal[nzone][j],rfield.flux[0][j]));
03833 }
03834 }
03835 }
03836
03837 else if( strcmp(save.chSave[ipPun],"SOUD") == 0 )
03838 {
03839
03840
03841 j = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon + 2;
03842 fprintf( save.ipPnunit[ipPun],
03843 "%.4e\t%.4e\t%.4e\t%.4e\n",
03844 opac.TauAbsFace[j-1],
03845 rfield.ConEmitLocal[nzone][j-1]/rfield.widflx[j-1]/MAX2(1e-35,opac.opacity_abs[j-1]),
03846 rfield.otscon[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1],
03847 rfield.otscon[iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]/opac.opacity_abs[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1] );
03848 }
03849
03850
03851 else if( strcmp(save.chSave[ipPun],"SPEC") == 0 )
03852 {
03853 SaveSpecial(save.ipPnunit[ipPun],chTime);
03854 }
03855
03856
03857 else if( strcmp(save.chSave[ipPun],"SPCS") == 0 )
03858 {
03859 if( ( strcmp(chTime,"LAST") != 0 && strcmp(save.chSaveArgs[ipPun],"COLU") != 0 ) ||
03860 ( strcmp(chTime,"LAST") == 0 && strcmp(save.chSaveArgs[ipPun],"COLU") == 0 ) )
03861 SaveSpecies(save.ipPnunit[ipPun] , ipPun);
03862 }
03863
03864 else if( strcmp(save.chSave[ipPun],"TEMP") == 0 )
03865 {
03866 static double deriv_old=-1;
03867 double deriv=-1. , deriv_sec;
03868
03869 fprintf( save.ipPnunit[ipPun], "%.5e\t%.4e\t%.2e",
03870 radius.depth_mid_zone,
03871 phycon.te,
03872 thermal.dCooldT );
03873
03874 if( nzone >1 )
03875 {
03876 deriv = (phycon.te - struc.testr[nzone-2])/ radius.drad;
03877 fprintf( save.ipPnunit[ipPun], "\t%.2e", deriv );
03878
03879 if( nzone > 2 )
03880 {
03881 deriv_sec = (deriv-deriv_old)/ radius.drad;
03882 fprintf( save.ipPnunit[ipPun], "\t%.2e",
03883 deriv_sec );
03884 }
03885 deriv_old = deriv;
03886 }
03887 fprintf( save.ipPnunit[ipPun], "\n");
03888 }
03889
03890
03891 else if( strcmp(save.chSave[ipPun],"TIMD") == 0 )
03892 {
03893 if( strcmp(chTime,"LAST") == 0 )
03894 DynaPunchTimeDep( save.ipPnunit[ipPun] , "END" );
03895 }
03896
03897
03898 else if( strcmp(save.chSave[ipPun],"XTIM") == 0 )
03899 {
03900 static double ElapsedTime , ZoneTime;
03901 if( nzone<=1 )
03902 {
03903 ElapsedTime = cdExecTime();
03904 ZoneTime = 0.;
03905 }
03906 else
03907 {
03908 double t = cdExecTime();
03909 ZoneTime = t - ElapsedTime;
03910 ElapsedTime = t;
03911 }
03912
03913
03914 fprintf( save.ipPnunit[ipPun], " %ld\t%.3f\t%.2f\n",
03915 nzone, ZoneTime , ElapsedTime );
03916 }
03917
03918 else if( strcmp(save.chSave[ipPun],"TPRE") == 0 )
03919 {
03920
03921 fprintf( save.ipPnunit[ipPun], "%5ld %11.4e %11.4e %11.4e %g\n",
03922 nzone, phycon.TeInit, phycon.TeProp, phycon.te,
03923 (phycon.TeProp- phycon.te)/phycon.te );
03924 }
03925
03926 else if( strcmp(save.chSave[ipPun],"WIND") == 0 )
03927 {
03928
03929
03930
03931 if( (save.punarg[ipPun][0] == 0 && strcmp(chTime,"LAST") == 0)
03932 ||
03933
03934 (save.punarg[ipPun][0] == 1 && strcmp(chTime,"LAST") != 0 ) )
03935 {
03936 fprintf( save.ipPnunit[ipPun],
03937 "%.5e\t%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
03938 radius.Radius_mid_zone,
03939 radius.depth_mid_zone,
03940 wind.windv,
03941 wind.AccelTotalOutward,
03942 wind.AccelLine,
03943 wind.AccelCont ,
03944 wind.fmul ,
03945 wind.AccelGravity );
03946 }
03947 }
03948
03949 else if( strcmp(save.chSave[ipPun],"XATT") == 0 )
03950 {
03951
03952 ASSERT( grid.lgOutputTypeOn[2] );
03953
03954 if( strcmp(chTime,"LAST") == 0 )
03955 {
03956 if( grid.lgGrid )
03957 saveFITSfile( save.ipPnunit[ipPun], 2 );
03958 else
03959 {
03960 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
03961 cdEXIT(EXIT_FAILURE);
03962 }
03963 }
03964 }
03965 else if( strcmp(save.chSave[ipPun],"XRFI") == 0 )
03966 {
03967
03968 ASSERT( grid.lgOutputTypeOn[3] );
03969
03970 if( strcmp(chTime,"LAST") == 0 )
03971 {
03972 if( grid.lgGrid )
03973 saveFITSfile( save.ipPnunit[ipPun], 3 );
03974 else
03975 {
03976 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
03977 cdEXIT(EXIT_FAILURE);
03978 }
03979 }
03980 }
03981 else if( strcmp(save.chSave[ipPun],"XINC") == 0 )
03982 {
03983
03984 ASSERT( grid.lgOutputTypeOn[1] );
03985
03986 if( strcmp(chTime,"LAST") == 0 )
03987 {
03988 if( grid.lgGrid )
03989 saveFITSfile( save.ipPnunit[ipPun], 1 );
03990 else
03991 {
03992 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
03993 cdEXIT(EXIT_FAILURE);
03994 }
03995 }
03996 }
03997 else if( strcmp(save.chSave[ipPun],"XDFR") == 0 )
03998 {
03999
04000 ASSERT( grid.lgOutputTypeOn[5] );
04001
04002 if( strcmp(chTime,"LAST") == 0 )
04003 {
04004 if( grid.lgGrid )
04005 saveFITSfile( save.ipPnunit[ipPun], 5 );
04006 else
04007 {
04008 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
04009 cdEXIT(EXIT_FAILURE);
04010 }
04011 }
04012 }
04013 else if( strcmp(save.chSave[ipPun],"XDFO") == 0 )
04014 {
04015
04016 ASSERT( grid.lgOutputTypeOn[4] );
04017
04018 if( strcmp(chTime,"LAST") == 0 )
04019 {
04020 if( grid.lgGrid )
04021 saveFITSfile( save.ipPnunit[ipPun], 4 );
04022 else
04023 {
04024 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
04025 cdEXIT(EXIT_FAILURE);
04026 }
04027 }
04028 }
04029 else if( strcmp(save.chSave[ipPun],"XLNR") == 0 )
04030 {
04031
04032 ASSERT( grid.lgOutputTypeOn[7] );
04033
04034 if( strcmp(chTime,"LAST") == 0 )
04035 {
04036 if( grid.lgGrid )
04037 saveFITSfile( save.ipPnunit[ipPun], 7 );
04038 else
04039 {
04040 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
04041 cdEXIT(EXIT_FAILURE);
04042 }
04043 }
04044 }
04045 else if( strcmp(save.chSave[ipPun],"XLNO") == 0 )
04046 {
04047
04048 ASSERT( grid.lgOutputTypeOn[6] );
04049
04050 if( strcmp(chTime,"LAST") == 0 )
04051 {
04052 if( grid.lgGrid )
04053 saveFITSfile( save.ipPnunit[ipPun], 6 );
04054 else
04055 {
04056 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
04057 cdEXIT(EXIT_FAILURE);
04058 }
04059 }
04060 }
04061 else if( strcmp(save.chSave[ipPun],"XREF") == 0 )
04062 {
04063
04064 ASSERT( grid.lgOutputTypeOn[9] );
04065
04066 if( strcmp(chTime,"LAST") == 0 )
04067 {
04068 if( grid.lgGrid )
04069 saveFITSfile( save.ipPnunit[ipPun], 9 );
04070 else
04071 {
04072 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
04073 cdEXIT(EXIT_FAILURE);
04074 }
04075 }
04076 }
04077 else if( strcmp(save.chSave[ipPun],"XTOT") == 0 )
04078 {
04079
04080 ASSERT( grid.lgOutputTypeOn[0] );
04081
04082 if( strcmp(chTime,"LAST") == 0 )
04083 {
04084 if( grid.lgGrid )
04085 saveFITSfile( save.ipPnunit[ipPun], 0 );
04086 else
04087 {
04088 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
04089 cdEXIT(EXIT_FAILURE);
04090 }
04091 }
04092 }
04093 else if( strcmp(save.chSave[ipPun],"XTRN") == 0 )
04094 {
04095
04096 ASSERT( grid.lgOutputTypeOn[8] );
04097
04098 if( strcmp(chTime,"LAST") == 0 )
04099 {
04100 if( grid.lgGrid )
04101 saveFITSfile( save.ipPnunit[ipPun], 8 );
04102 else
04103 {
04104 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
04105 cdEXIT(EXIT_FAILURE);
04106 }
04107 }
04108 }
04109 else if( strcmp(save.chSave[ipPun],"XSPM") == 0 )
04110 {
04111
04112 ASSERT( grid.lgOutputTypeOn[10] );
04113
04114 if( strcmp(chTime,"LAST") == 0 )
04115 {
04116 if( grid.lgGrid )
04117 saveFITSfile( save.ipPnunit[ipPun], 10 );
04118 else
04119 {
04120 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
04121 cdEXIT(EXIT_FAILURE);
04122 }
04123 }
04124 }
04125
04126
04127
04128
04129
04130 else if( save.lgRealSave[ipPun] && lgNoHitFirstBranch )
04131 {
04132
04133 fprintf( ioQQQ, " PROBLEM DISASTER SaveDo does not recognize flag %4.4s set by ParseSave. This is impossible.\n",
04134 save.chSave[ipPun] );
04135 TotalInsanity();
04136 }
04137
04138
04139
04140
04141
04142
04143
04144
04145
04146 if( strcmp(chTime,"LAST") == 0 &&
04147 !(iterations.lgLastIt && !grid.lgGrid ) &&
04148 save.lgHashEndIter[ipPun] &&
04149 save.lg_separate_iterations[ipPun] &&
04150 !save.lgFITS[ipPun] )
04151 {
04152 if( dynamics.lgTimeDependentStatic && strcmp( save.chHashString , "TIME_DEP" )==0 )
04153 {
04154 fprintf( save.ipPnunit[ipPun], "\"time=%f\n",
04155 dynamics.time_elapsed );
04156 }
04157 else
04158 {
04159 fprintf( save.ipPnunit[ipPun], "%s",
04160 save.chHashString );
04161 if( grid.lgGrid && ( iterations.lgLastIt || lgAbort ) )
04162 fprintf( save.ipPnunit[ipPun], " GRID_DELIMIT -- grid%09ld",
04163 optimize.nOptimiz );
04164 fprintf( save.ipPnunit[ipPun], "\n" );
04165 }
04166 }
04167 if( save.lgFLUSH )
04168 fflush( save.ipPnunit[ipPun] );
04169 }
04170 }
04171 return;
04172 }
04173
04174
04175 STATIC void SaveLineIntensity(FILE * ioPUN, long int ipPun , realnum Threshold )
04176 {
04177 long int i;
04178
04179 DEBUG_ENTRY( "SaveLineIntensity()" );
04180
04181
04182
04183
04184 fprintf( ioPUN, "**********************************************************************************************************************************\n" );
04185 input.echo(ioPUN);
04186
04187
04188 cdWarnings( ioPUN);
04189 cdCautions( ioPUN);
04190 fprintf( ioPUN, "zone=%5ld\n", nzone );
04191 fprintf( ioPUN, "**********************************************************************************************************************************\n" );
04192 fprintf( ioPUN, "begin emission lines\n" );
04193
04194
04195 SaveResults1Line(ioPUN," ",0,0.,"Start");
04196
04197
04198 bool lgEmergent = false;
04199 if( save.punarg[ipPun][0] > 0 )
04200 lgEmergent = true;
04201
04202 for( i=0; i < LineSave.nsum; i++ )
04203 {
04204
04205
04206 if( LineSv[i].SumLine[lgEmergent] > Threshold )
04207 {
04208 SaveResults1Line(ioPUN,(char*)LineSv[i].chALab,LineSv[i].wavelength,
04209 LineSv[i].SumLine[save.lgEmergent[ipPun]], "Line ");
04210 }
04211 }
04212
04213 SaveResults1Line(ioPUN," ",0,0.,"Flush");
04214
04215 fprintf( ioPUN, " \n" );
04216 fprintf( ioPUN, "**********************************************************************************************************************************\n" );
04217
04218 return;
04219 }
04220
04221
04222 static bool lgPopsFirstCall , lgSaveOpticalDepths;
04223
04224
04225 STATIC void SaveLineStuff(
04226 FILE * ioPUN,
04227 const char *chJob ,
04228 realnum xLimit )
04229 {
04230
04231 long int nelem , ipLo , ipHi , i , ipISO;
04232 long index = 0;
04233 static bool lgFirst=true;
04234
04235 DEBUG_ENTRY( "SaveLineStuff()" );
04236
04237
04238 if( strcmp( &*chJob , "optical" ) == 0 )
04239 {
04240
04241 lgSaveOpticalDepths = true;
04242 lgPopsFirstCall = false;
04243 }
04244 else if( strcmp( &*chJob , "populat" ) == 0 )
04245 {
04246 lgSaveOpticalDepths = false;
04247
04248 if( lgFirst )
04249 {
04250 lgPopsFirstCall = true;
04251 fprintf(ioPUN,"index\tAn.ion\tgLo\tgUp\tE(wn)\tgf\n");
04252 lgFirst = false;
04253 }
04254 else
04255 {
04256 lgPopsFirstCall = false;
04257 }
04258 }
04259 else
04260 {
04261 fprintf( ioQQQ, " insane job in SaveLineStuff =%s\n",
04262 &*chJob );
04263 cdEXIT(EXIT_FAILURE);
04264 }
04265
04266
04267
04268
04269 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
04270 {
04271 for( nelem=ipISO; nelem < LIMELM; nelem++ )
04272 {
04273 if( dense.lgElmtOn[nelem] )
04274 {
04275
04276 for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
04277 {
04278 for( ipLo=0; ipLo <ipHi; ipLo++ )
04279 {
04280 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
04281 continue;
04282
04283 ++index;
04284 Save1Line( iso_sp[ipISO][nelem].trans(ipHi,ipLo), ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[nelem]) );
04285 }
04286 }
04287
04288
04289
04290 if( lgSaveOpticalDepths )
04291 {
04292
04293
04294
04295
04296
04297 for( ipHi=iso_sp[ipISO][nelem].st[iso_sp[ipISO][nelem].numLevels_local-1].n()+1; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )
04298 {
04299 ++index;
04300 Save1Line( ExtraLymanLines[ipISO][nelem][ipExtraLymanLines[ipISO][nelem][ipHi]], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[nelem]) );
04301 }
04302 }
04303 }
04304 }
04305 }
04306
04307
04308 for( i=1; i < nLevel1; i++ )
04309 {
04310 ++index;
04311 Save1Line( TauLines[i], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[(*TauLines[i].Hi()).nelem()-1]) );
04312 }
04313
04314 for( i=0; i < nWindLine; i++ )
04315 {
04316 if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
04317 {
04318 ++index;
04319 Save1Line( TauLine2[i], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[(*TauLine2[i].Hi()).nelem()-1]) );
04320 }
04321 }
04322
04323 for( i=0; i < nUTA; i++ )
04324 {
04325 ++index;
04326 Save1Line( UTALines[i], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[(*UTALines[i].Hi()).nelem()-1]) );
04327 }
04328
04329
04330
04331 FeIIPunchLineStuff( ioPUN , xLimit , index);
04332
04333
04334 h2.H2_PunchLineStuff( ioPUN , xLimit , index);
04335
04336
04337 fprintf( ioPUN , "%s\n",save.chHashString );
04338 return;
04339 }
04340
04341
04342 void Save1Line( const TransitionProxy& t , FILE * ioPUN , realnum xLimit , long index, realnum DopplerWidth )
04343 {
04344
04345 if( lgSaveOpticalDepths )
04346 {
04347
04348 if( t.Emis().TauIn() >= xLimit )
04349 {
04350
04351 fprintf( ioPUN , "%2.2s%2.2s\t",
04352 elementnames.chElementSym[(*t.Hi()).nelem()-1] ,
04353 elementnames.chIonStage[(*t.Hi()).IonStg()-1] );
04354
04355
04356
04357 if( strcmp( save.chConPunEnr[save.ipConPun], "labl" )== 0 )
04358 {
04359 prt_wl( ioPUN , t.WLAng() );
04360 }
04361 else
04362 {
04363
04364 fprintf( ioPUN , "%.7e", AnuUnit((realnum)(t.EnergyRyd())) );
04365 }
04366
04367 fprintf( ioPUN , "\t%.3f", t.Emis().TauIn() );
04368
04369 fprintf(ioPUN, "\t%.3e",
04370 t.Emis().dampXvel() / DopplerWidth );
04371 fprintf(ioPUN, "\n");
04372 }
04373 }
04374 else if( lgPopsFirstCall )
04375 {
04376 char chLbl[11];
04377
04378 strcpy( chLbl, chLineLbl(t) );
04379 fprintf(ioPUN, "%li\t%s" , index , chLbl );
04380
04381 fprintf(ioPUN, "\t%.0f\t%.0f",
04382 (*t.Lo()).g() ,(*t.Hi()).g());
04383
04384 fprintf(ioPUN, "\t%.2f\t%.3e",
04385 t.EnergyWN() ,t.Emis().gf());
04386 fprintf(ioPUN, "\n");
04387 }
04388 else
04389 {
04390
04391 if( (*t.Hi()).Pop() > xLimit )
04392 {
04393
04394
04395
04396
04397
04398 fprintf(ioPUN,"%li\t%.2e\t%.2e\n", index, (*t.Lo()).Pop(), (*t.Hi()).Pop() );
04399 }
04400 }
04401 }
04402
04403
04404 STATIC void SaveNewContinuum(FILE * ioPUN )
04405 {
04406 long int ipLo, ipHi,
04407 j ,
04408 ncells;
04409
04410 double wllo, wlhi;
04411
04412 double *energy,
04413 *cont_incid,
04414 *cont_atten,
04415 *diffuse_in,
04416 *diffuse_out,
04417 *emis_lines_out,
04418 *emis_lines_in;
04419
04420
04421 wllo = rfield.anu[0];
04422
04423 wlhi = rfield.anu[rfield.nflux-1];
04424
04425 ipLo = ipoint(wllo)-1;
04426 ipHi = ipoint(wlhi)-1;
04427 ncells = ipHi - ipLo + 1;
04428
04429
04430 energy = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
04431 cont_incid = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
04432 cont_atten = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
04433 diffuse_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
04434 diffuse_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
04435 emis_lines_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
04436 emis_lines_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
04437
04438
04439
04440
04441 for( j=0; j<ncells; ++j )
04442 {
04443 energy[j] = rfield.AnuOrg[j+ipLo] - rfield.widflx[j+ipLo]/2.;
04444 }
04445
04446 fixit();
04447
04448
04449
04450
04451 #if 1
04452
04453
04454 cdSPEC( 1 , ncells , cont_incid );
04455
04456 cdSPEC( 2 , ncells , cont_atten );
04457
04458 cdSPEC( 5 , ncells , diffuse_in );
04459
04460 cdSPEC( 4 , ncells , diffuse_out );
04462
04463 cdSPEC( 6 , ncells , emis_lines_out );
04464
04465 cdSPEC( 7 , ncells , emis_lines_in );
04466 #else
04467 cdSPEC2( 1, rfield.nflux, 0, rfield.nflux - 1, cont_incid );
04468
04469 cdSPEC2( 2, rfield.nflux, 0, rfield.nflux - 1, cont_atten );
04470
04471 cdSPEC2( 5, rfield.nflux, 0, rfield.nflux - 1, diffuse_in );
04472
04473 cdSPEC2( 4, rfield.nflux, 0, rfield.nflux - 1, diffuse_out );
04474
04475 cdSPEC2( 6, rfield.nflux, 0, rfield.nflux - 1, emis_lines_out );
04476
04477 cdSPEC2( 7, rfield.nflux, 0, rfield.nflux - 1, emis_lines_in );
04478 #endif
04479
04480
04481 for( j=0; j<ncells-1; ++j )
04482 {
04483
04484 fprintf( ioPUN,"%.5e\t", AnuUnit((realnum)(energy[j]+rfield.widflx[j+ipLo]/2.) ) );
04485 fprintf( ioPUN,"%.3e\t", cont_incid[j] );
04486 fprintf( ioPUN,"%.3e\t", cont_atten[j] );
04487 fprintf( ioPUN,"%.3e\t", diffuse_in[j]+diffuse_out[j] );
04488 fprintf( ioPUN,"%.3e",
04489 emis_lines_out[j]+emis_lines_in[j] );
04490 fprintf( ioPUN,"\n" );
04491 }
04492
04493 free(energy);
04494 free(cont_incid);
04495 free(diffuse_in);
04496 free(diffuse_out);
04497 free(cont_atten);
04498 free(emis_lines_out);
04499 free(emis_lines_in);
04500
04501
04502 }
04503
04504
04505 STATIC void AGN_Hemis(FILE *ioPUN )
04506 {
04507 const int NTE = 4;
04508 realnum te[NTE] = {5000., 10000., 15000., 20000.};
04509 realnum *agn_continuum[NTE];
04510 double TempSave = phycon.te;
04511 long i , j;
04512
04513 DEBUG_ENTRY( "AGN_Hemis()" );
04514
04515
04516
04517 for( i=0;i<NTE; ++i)
04518 {
04519 agn_continuum[i] = (realnum *)MALLOC((unsigned)rfield.nflux*sizeof(realnum) );
04520
04521
04522
04523 TempChange(te[i] , true);
04524
04525 ConvPresTempEdenIoniz();
04526
04527
04528 RT_diffuse();
04529 for(j=0;j<rfield.nflux; ++j )
04530 {
04531 agn_continuum[i][j] = rfield.ConEmitLocal[nzone][j]/(realnum)dense.eden/
04532 (dense.xIonDense[ipHYDROGEN][1] + dense.xIonDense[ipHELIUM][1] + dense.xIonDense[ipHELIUM][2] );
04533 }
04534 }
04535
04536
04537 fprintf(ioPUN,"wl");
04538 for( i=0;i<NTE; ++i)
04539 {
04540 fprintf(ioPUN,"\tT=%.0f",te[i]);
04541 }
04542 fprintf( ioPUN , "\tcont\n");
04543
04544
04545 for(j=0;j<rfield.nflux; ++j )
04546 {
04547 fprintf( ioPUN , "%12.5e",
04548 AnuUnit(rfield.AnuOrg[j]) );
04549
04550 for( i=0;i<NTE; ++i)
04551 {
04552 fprintf(ioPUN,"\t%.3e",agn_continuum[i][j]*rfield.anu2[j]*EN1RYD/rfield.widflx[j]);
04553 }
04554
04555 fprintf( ioPUN , "\t%s\n" , rfield.chContLabel[j]);
04556 }
04557
04558
04559 for( i=0;i<NTE; ++i)
04560 {
04561 free( agn_continuum[i] );
04562 }
04563
04564
04565
04566 TempChange(TempSave , true);
04567
04568 fprintf( ioQQQ, "AGN_Hemis - result of save AGN3 hemis - I have left the code in a disturbed state, and will now exit.\n");
04569 cdEXIT(EXIT_FAILURE);
04570 }
04571
04572
04573
04574 STATIC void SaveResults(FILE* ioPUN)
04575 {
04576 long int i , nelem , ion;
04577
04578 DEBUG_ENTRY( "SaveResults()" );
04579
04580
04581
04582
04583 fprintf( ioPUN, "**********************************************************************************************************************************\n" );
04584 input.echo(ioPUN);
04585
04586
04587 cdWarnings(ioPUN);
04588 cdCautions(ioPUN);
04589 fprintf( ioPUN, "**********************************************************************************************************************************\n" );
04590
04591 fprintf( ioPUN, "C*OPTICAL DEPTHS ELECTRON=%10.3e\n", opac.telec );
04592
04593 fprintf( ioPUN, "BEGIN EMISSION LINES\n" );
04594 SaveResults1Line(ioPUN," ",0,0.,"Start");
04595
04596 for( i=0; i < LineSave.nsum; i++ )
04597 {
04598 if( LineSv[i].SumLine[0] > 0. )
04599 {
04600 SaveResults1Line(ioPUN,(char*)LineSv[i].chALab,LineSv[i].wavelength,
04601 LineSv[i].SumLine[0],"Line ");
04602 }
04603 }
04604
04605 SaveResults1Line(ioPUN," ",0,0.,"Flush");
04606
04607 fprintf( ioPUN, " \n" );
04608
04609 fprintf( ioPUN, "BEGIN COLUMN DENSITIES\n" );
04610
04611
04612
04613
04614 ASSERT( LIMELM == 30 );
04615
04616
04617 for( nelem=0; nelem<LIMELM; nelem++ )
04618 {
04619 for(ion=0; ion < nelem+1; ion++)
04620 {
04621 fprintf( ioPUN, " %10.3e", mean.xIonMean[0][nelem][ion][0] );
04622
04623 if( nelem==9|| nelem==19 || nelem==29 )
04624 {
04625 fprintf( ioPUN, "\n" );
04626 }
04627 }
04628 }
04629
04630 fprintf( ioPUN, "END OF RESULTS\n" );
04631 fprintf( ioPUN, "**********************************************************************************************************************************\n" );
04632 return;
04633 }
04634
04635 static const int LINEWIDTH = 6;
04636
04637
04638
04639 STATIC void SaveResults1Line(
04640 FILE * ioPUN,
04641
04642 const char *chLab,
04643 realnum wl,
04644 double xInten,
04645
04646 const char *chFunction)
04647 {
04648
04649 long int i;
04650 static realnum wavelength[LINEWIDTH];
04651 static long ipLine;
04652 static double xIntenSave[LINEWIDTH];
04653 static char chLabSave[LINEWIDTH][5];
04654
04655 DEBUG_ENTRY( "SaveResults1Line()" );
04656
04657
04658
04659 if( strcmp(chFunction,"Start") == 0 )
04660 {
04661 ipLine = 0;
04662 }
04663 else if( strcmp(chFunction,"Line ") == 0 )
04664 {
04665
04666 wavelength[ipLine] = wl;
04667 strcpy( chLabSave[ipLine], chLab );
04668 xIntenSave[ipLine] = xInten;
04669
04670
04671
04672 ++ipLine;
04673
04674
04675 if( ( strcmp(save.chPunRltType,"column") == 0 ) || ipLine == LINEWIDTH )
04676 {
04677
04678 for( i=0; i < ipLine; i++ )
04679 {
04680 fprintf( ioPUN, " %4.4s ", chLabSave[i] );
04681 prt_wl( ioPUN, wavelength[i] );
04682 fprintf( ioPUN,"\t%.3e", xIntenSave[i]);
04683
04684
04685 if( strcmp(save.chPunRltType,"column") == 0 )
04686 fprintf( ioPUN, "\n" );
04687 }
04688
04689 if( strcmp(save.chPunRltType,"array ") == 0 )
04690 fprintf( ioPUN, " \n" );
04691 ipLine = 0;
04692 }
04693 }
04694 else if( strcmp(chFunction,"Flush") == 0 )
04695 {
04696 if( ipLine > 0 )
04697 {
04698
04699
04700
04701
04702
04703 for( i=0; i < ipLine; i++ )
04704 {
04705 fprintf( ioPUN, " %4.4s", chLabSave[i] );
04706 prt_wl( ioPUN, wavelength[i] );
04707 fprintf( ioPUN,"\t%.3e", xIntenSave[i]);
04708
04709
04710 if( strcmp(save.chPunRltType,"column") == 0 )
04711 fprintf( ioPUN, "\n" );
04712 }
04713 if( strcmp(save.chPunRltType,"array ") == 0 )
04714 fprintf( ioPUN, " \n" );
04715 }
04716 }
04717 else
04718 {
04719 fprintf( ioQQQ, " SaveResults1Line called with insane request=%5.5s\n",
04720 chFunction );
04721 cdEXIT(EXIT_FAILURE);
04722 }
04723 return;
04724 }
04725
04726 static const int NENR_GAUNT = 37;
04727 static const int NTE_GAUNT = 21;
04728
04729
04730 STATIC void SaveGaunts(FILE* ioPUN)
04731 {
04732 long int i,
04733 charge,
04734 ite,
04735 j;
04736
04737 realnum ener[NENR_GAUNT],
04738 ste[NTE_GAUNT],
04739 z;
04740 double tempsave;
04741 double g[NENR_GAUNT][NTE_GAUNT], gfac;
04742
04743
04744 DEBUG_ENTRY( "SaveGaunts()" );
04745
04746
04747
04748 tempsave = phycon.te;
04749
04750 for( i=0; i < NTE_GAUNT; i++ )
04751 {
04752 ste[i] = 0.5f*i;
04753 }
04754
04755 for( i=0; i < NENR_GAUNT; i++ )
04756 {
04757 ener[i] = 0.5f*i - 8.f;
04758 }
04759
04760 for( charge = 1; charge<LIMELM; charge++ )
04761 {
04762
04763 z = (realnum)log10((double)charge);
04764
04765
04766 for( ite=0; ite < NTE_GAUNT; ite++ )
04767 {
04768 phycon.alogte = ste[ite];
04769 phycon.te = pow(10.,phycon.alogte);
04770
04771 for( j=0; j < NENR_GAUNT; j++ )
04772 {
04773 gfac = cont_gaunt_calc( phycon.te, (double)charge, pow(10.,(double)ener[j]) );
04774
04775
04776 g[j][ite] = gfac;
04777 }
04778 }
04779
04780
04781 fprintf( ioPUN, " " );
04782 for( i=1; i <= NTE_GAUNT; i++ )
04783 {
04784 fprintf( ioPUN, "\t%6.1f", ste[i-1] );
04785 }
04786 fprintf( ioPUN, "\n" );
04787
04788 for( j=0; j < NENR_GAUNT; j++ )
04789 {
04790 fprintf( ioPUN, "\t%6.2f", ener[j] );
04791 for( ite=0; ite < NTE_GAUNT; ite++ )
04792 {
04793 fprintf( ioPUN, "\t%6.2f", log10(g[j][ite]) );
04794 }
04795 fprintf( ioPUN, "\n" );
04796 }
04797
04798 fprintf( ioPUN, " " );
04799 for( i=0; i < NTE_GAUNT; i++ )
04800 {
04801 fprintf( ioPUN, "\t%6.1f", ste[i] );
04802 }
04803 fprintf( ioPUN, "\n" );
04804
04805
04806 fprintf( ioPUN, " " );
04807 for( i=0; i < NTE_GAUNT; i++ )
04808 {
04809 fprintf( ioPUN, "\t%6.1f", ste[i] );
04810 }
04811 fprintf( ioPUN, "\n" );
04812
04813 for( i=0; i < NTE_GAUNT; i++ )
04814 {
04815 for( j=0; j < NENR_GAUNT; j++ )
04816 {
04817 fprintf( ioPUN, "\t%6.2f\t%6.2f\t%6.2e\n", ste[i], ener[j],
04818 log10(g[j][i]) );
04819 }
04820 }
04821
04822 fprintf( ioPUN, "Below is log(u), log(gamma**2), gff\n" );
04823
04824 for( i=0; i < NTE_GAUNT; i++ )
04825 {
04826 for( j=0; j < NENR_GAUNT; j++ )
04827 {
04828 fprintf( ioPUN, "\t%6.2f\t%6.2f\t%6.2e\n", 2.*z + log10(TE1RYD) - ste[i] , log10(TE1RYD)+ener[j]-ste[i],
04829 log10(g[j][i]) );
04830 }
04831 }
04832 fprintf( ioPUN, "end of charge = %li\n", charge);
04833 fprintf( ioPUN, "****************************\n");
04834 }
04835
04836 phycon.te = tempsave;
04837 phycon.alogte = log10( phycon.te );
04838 return;
04839 }
04840
04841 void SaveGrid(FILE* pnunit, exit_type status)
04842 {
04843 if( pnunit == NULL )
04844 return;
04845
04846 if( optimize.nOptimiz == 0 )
04847 {
04848
04849 fprintf( pnunit, "#Index\tFailure?\tWarnings?\tExit code\t#rank\t#seq" );
04850
04851 for( int i=0; i < grid.nintparm; i++ )
04852 {
04853 char chStr[10];
04854 strncpy( chStr, optimize.chVarFmt[i], 9 );
04855
04856 chStr[9] = '\0';
04857 fprintf( pnunit, "\t%s", chStr );
04858 }
04859 fprintf( pnunit, "\tgrid parameter string\n" );
04860 }
04861
04862 bool lgNoFailure = ( status == ES_SUCCESS || status == ES_WARNINGS );
04863 fprintf( pnunit, "%9.9ld\t%c\t%c\t%20s\t%ld\t%ld",
04864 optimize.nOptimiz,
04865 TorF(!lgNoFailure),
04866 TorF(warnings.lgWarngs),
04867 cpu.i().chExitStatus(status).c_str(),
04868 cpu.i().nRANK(),
04869 grid.seqNum );
04870
04871 char chGridParam[INPUT_LINE_LENGTH];
04872 char chStringHold[100];
04873 sprintf( chStringHold, "%f", grid.interpParameters[optimize.nOptimiz][0] );
04874 strcpy( chGridParam, chStringHold );
04875 for( int j=0; j < grid.nintparm; j++ )
04876 {
04877 if( j > 0 )
04878 {
04879 sprintf( chStringHold, ", %f", grid.interpParameters[optimize.nOptimiz][j] );
04880 strcat( chGridParam, chStringHold );
04881 }
04882 fprintf( pnunit, "\t%f", grid.interpParameters[optimize.nOptimiz][j] );
04883 }
04884 fprintf( pnunit, "\t%s\n", chGridParam );
04885 }