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