/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and * others. For conditions of distribution and use see copyright notice in license.txt */ /*SaveDo produce save output during calculation, * chTime is 'MIDL' during calculation, 'LAST' at the end */ /*SaveNewContinuum produce the 'save new continuum' output */ /*SaveLineStuff save optical depths or source functions for all transferred lines */ /*Save1Line called by SaveLineStuff to produce output for one line */ /*SaveNewContinuum produce the 'save new continuum' output */ /*SaveLineIntensity produce the 'save lines intensity' output */ /* save h emission, for AGN3 chapter 4, routine is below */ /*SaveResults1Line do single line of output for the save results and save line intensity commands */ /* the number of emission lines across one line of printout */ /*SaveSpecial generate output for the save special command */ /*SaveResults save results from save results command */ /*SaveResults1Line do single line of output for the save results and save line intensity commands */ /*SaveGaunts called by save gaunts command to output gaunt factors */ /*SaveFeII_cont return appropriate FeII banded continuum, inward, outward, or total */ /*FindStrongestLineLabels find strongest lines contributing to point in continuum energy mesh, output in some save commands */ #include "cddefines.h" #include "cddrive.h" #include "physconst.h" #include "mean.h" #include "taulines.h" #include "struc.h" #include "iso.h" #include "mole.h" #include "hyperfine.h" #include "rt.h" #include "lines_service.h" #include "doppvel.h" #include "dense.h" #include "h2.h" #include "magnetic.h" #include "hydrogenic.h" #include "secondaries.h" #include "grainvar.h" #include "lines.h" #include "dynamics.h" #include "colden.h" #include "continuum.h" #include "ionbal.h" #include "yield.h" #include "prt.h" #include "iterations.h" #include "heavy.h" #include "conv.h" #include "geometry.h" #include "called.h" #include "helike.h" #include "opacity.h" #include "rfield.h" #include "phycon.h" #include "timesc.h" #include "radius.h" #include "atomfeii.h" #include "monitor_results.h" #include "thermal.h" #include "wind.h" #include "hmi.h" #include "pressure.h" #include "elementnames.h" #include "ipoint.h" #include "gammas.h" #include "atmdat.h" #include "hcmap.h" #include "input.h" #include "save.h" #include "optimize.h" #include "warnings.h" #include "grid.h" #include "mole_priv.h" #include "atomfeii.h" // this needs C linkage since it is passed as a pointer to a C library routine extern "C" int wavelength_compare (const void * a, const void * b) { LinSv a1 = *(LinSv*)a; LinSv b1 = *(LinSv*)b; /* comparison is b-a so we get inverse wavelength order (increasing energy order) */ if( b1.wavelength > a1.wavelength ) return 1; else if( b1.wavelength < a1.wavelength ) return -1; else return 0; } // find strongest lines contributing to point in continuum energy mesh, output in some save commands STATIC void FindStrongestLineLabels( void ) { long low_index=0; long high_index=0; long j_min = 0; double MaxFlux = 0.; long ipMaxFlux = 0; long j = 0; ASSERT( LineSave.ipass==1 ); // copy and sort memcpy(LineSvSortWL, LineSv, (unsigned)LineSave.nsum*sizeof( LinSv )); qsort((void *)LineSvSortWL,(size_t)LineSave.nsum,sizeof( LinSv ),wavelength_compare); while( rfield.anu[j_min]+0.5*rfield.widflx[j_min] < RYDLAM/LineSvSortWL[0].wavelength) j_min++; for( j=0; j 0. ); ASSERT( LineSvSortWL[high_index].wavelength > 0. ); ASSERT( RYDLAM/LineSvSortWL[low_index].wavelength >= rfield.anu[j]-0.5*rfield.widflx[j] ); ASSERT( RYDLAM/LineSvSortWL[high_index].wavelength <= rfield.anu[j]+0.5*rfield.widflx[j] ); MaxFlux = 0.; ipMaxFlux = 0; for( long k = low_index; k <= high_index; k++ ) { if( strcmp( LineSvSortWL[k].chALab, "Coll" )==0 || strcmp( LineSvSortWL[k].chALab, "Heat" )==0 || strcmp( LineSvSortWL[k].chALab, "Pump" )==0 || strcmp( LineSvSortWL[k].chALab, "pump" )==0 || strcmp( LineSvSortWL[k].chALab, "nInu" )==0 || strcmp( LineSvSortWL[k].chALab, "nFnu" )==0 || strcmp( LineSvSortWL[k].chALab, "InwT" )==0 || strcmp( LineSvSortWL[k].chALab, "InwC" )==0 || strcmp( LineSvSortWL[k].chALab, "Inwd" )==0 || strcmp( LineSvSortWL[k].chALab, "Ca A" )==0 || strcmp( LineSvSortWL[k].chALab, "Ca B" )==0 || strcmp( LineSvSortWL[k].chALab, "Pho+" )==0 || strcmp( LineSvSortWL[k].chALab, "Pcon" )==0 || strcmp( LineSvSortWL[k].chALab, "Q(H)" )==0 || strcmp( LineSvSortWL[k].chALab, "Unit" )==0 ) continue; if( LineSvSortWL[k].SumLine[0] > MaxFlux ) { MaxFlux = LineSvSortWL[k].SumLine[0]; ipMaxFlux = k; } } /* line label */ if( ipMaxFlux > 0 ) strcpy( rfield.chLineLabel[j], LineSvSortWL[ipMaxFlux].chALab ); } return; } # ifdef USE_NLTE7 //Run "Save NLTE" static void runNLTE(int ipPun) { long nelem = ipNEON; int nouter[11] = {0,2,3,3,3,4,3,5,5,6,0}; ASSERT(dense.gas_phase[nelem] != 0.); // Section 1 fprintf( save.ipPnunit[ipPun], "data Ferland University of Kentucky Cloudy 10 11-5-2011\n"); fprintf( save.ipPnunit[ipPun], "case NeXY6\n"); fprintf( save.ipPnunit[ipPun], "code Cloudy\n"); fprintf( save.ipPnunit[ipPun], "atom Ne 10\n"); fprintf( save.ipPnunit[ipPun], "calctime %11.4e %11.4e\n",4.5,30.); //Section 2 /**Calculations**/ //Avg Charge double avgchrg = 0.; double m2 = 0.; double m3 = 0.; //double m4 = 0.; for( int ion = 1; ion < nelem-1; ion++) { avgchrg += (dense.xIonDense[nelem][ion]/dense.gas_phase[nelem] * ion); } //Central Moments of Charge Distribution for( int ion = 1; ion < nelem-1; ion++) { m2 += dense.xIonDense[nelem][ion]/dense.gas_phase[nelem]*POW2((ion - avgchrg)); m3 += dense.xIonDense[nelem][ion]/dense.gas_phase[nelem]*POW3((ion - avgchrg)); //m4 += dense.xIonDense[nelem][ion]/dense.gas_phase[nelem]*POW4((ion - avgchrg)); } //Specific internal energy double Eint = 0.; double partfun = 0.; long count_nrglvls = 0; for( int ion = 1; ion < nelem-1; ion++) { // Find ipSpecies long ipSpec = 0; for( int i = 0; i < nSpecies; i++) { ipSpec = i; if(dBaseStates[i][0].nelem()-1 == nelem && dBaseStates[i][0].IonStg()-1 == ion) break; ipSpec = -1; } if( ipSpec != -1) {//This requires a loop over energy levels for( int lvl = 0; lvl < dBaseSpecies[ipSpec].numLevels_max; lvl++) { Eint += dBaseStates[ipSpec][lvl].Pop()*dBaseStates[ipSpec][lvl].energy().eV(); partfun += dBaseStates[ipSpec][lvl].g()*exp(-1*dBaseStates[ipSpec][lvl].energy().eV()/ phycon.te_eV); count_nrglvls += 1; } } else { printf("\nNot using :\t%li\t%i\n",nelem+1,ion+1); //cdEXIT(EXIT_FAILURE); } } Eint = Eint / dense.gas_phase[nelem]; fprintf( save.ipPnunit[ipPun],"\nsummary_quantities\n"); fprintf( save.ipPnunit[ipPun],"plasma %11.4e %11.4e\n",phycon.te_eV,dense.eden); fprintf( save.ipPnunit[ipPun],"time %11.4e\n",0.); fprintf( save.ipPnunit[ipPun],"zbar %11.4e\n",avgchrg); fprintf( save.ipPnunit[ipPun],"m2 %11.4e\n",m2); fprintf( save.ipPnunit[ipPun],"m3 %11.4e\n",m3); //fprintf( save.ipPnunit[ipPun],"m4 %11.3e\n",m4); fprintf( save.ipPnunit[ipPun],"eint %11.4e\n",0.); //%11.3e\n",Eint); fprintf( save.ipPnunit[ipPun],"deintdt %11.4e\n",0.); //%11.4e\n",thermal.dCooldT); fprintf( save.ipPnunit[ipPun],"pfn %11.4e\n",partfun); fprintf( save.ipPnunit[ipPun],"nmax_eff %i\n",0); fprintf( save.ipPnunit[ipPun],"ploss %11.4e %11.4e %11.4e %9.3e\n",0.,0.,0.,thermal.totcol); //Section 3 fprintf( save.ipPnunit[ipPun],"\nion_stages %li\n",nelem-2); //Photoionization Rate double PIR[nelem+1]; //autoionization Rate double autoion[nelem+1]; //This requires a loop over ion stages for( int ion = 1; ion < nelem-1; ion++) { double tempRecomb = 0.; double f_aauto = 0.; double f_aphoto = 0.; double f_acoll = 0.; double f_Scoll = 0.; double f_Sphoto = 0.; double f_auto = 0.; autoion[ion] = 0.; PIR[ion] = 0.; //This deals with the fact that there are no recombinations FROM atoms if( ion != 0 && ionbal.RateRecomTot[nelem][ion-1] != 0.) { tempRecomb = ionbal.RateRecomTot[nelem][ion-1]; f_aauto = dense.eden*ionbal.DR_Badnell_rate_coef[nelem][ion-1]/tempRecomb; f_aphoto = dense.eden*ionbal.RR_rate_coef_used[nelem][ion-1]/tempRecomb; f_acoll = dense.eden*ionbal.CotaRate[ion-1]/tempRecomb; } //Photoionization rate for( int i=0; i < Heavy.nsShells[nelem][ion]; i++) { //Add up the photoionization rates for all of the shells if( ion != nelem+1) { PIR[ion] += ionbal.PhotoRate_Shell[nelem][ion][i][0]; } } //Compute the fractional ionizations if( ion != nelem+1 && ionbal.RateIonizTot(nelem,ion) != 0.) { f_Scoll = ionbal.CollIonRate_Ground[nelem][ion][0]/ionbal.RateIonizTot(nelem,ion); f_Sphoto = PIR[ion]/ionbal.RateIonizTot(nelem,ion); f_auto = autoion[ion]/ionbal.RateIonizTot(nelem,ion); } //Autoionizations autoion[ion] = ionbal.UTA_ionize_rate[nelem][ion] + secondaries.csupra[nelem][ion]; if( ion > 0) { autoion[ion] += 0.0; //auger[ion-1]; } fprintf( save.ipPnunit[ipPun],"ion %li %11.4e %i\n" " %11.3e %11.3e %11.3e %11.3e\n" " %11.3e %11.3e %11.3e %11.3e\n\n", nelem+1-ion, dense.xIonDense[nelem][ion]/dense.gas_phase[nelem], nouter[ion], ionbal.RateIonizTot(nelem,ion), f_Scoll, f_Sphoto, f_auto, tempRecomb, f_acoll, f_aphoto, f_aauto); } //Section 4 fprintf( save.ipPnunit[ipPun],"\nenergy_levels %li\n",count_nrglvls); //Collisional Destruction Rate Bound-Free - Collisional Ionization double GcollBF = 0.0; //Photo Dest Rate Bound-Free - Photoionization double GphotoBF = 0.0; //Autoionization Dest Rate double GautoBF = 0.0; // Total Dest Rate double Gtotal = 0.0; //Collisional Creation Rate Bound-Free - Collisional Ionization double QcollBF = 0.0; //Photo Creation Rate Bound-Free - Photoionization double QphotoBF = 0.0; //Autoionization Creation Rate double QautoBF = 0.0; // Total Creation Rate double Qtotal = 0.0; //total pop of all levels double totallevelpop = 0.; for( int ion = 1; ion < nelem-1; ion++) { // Find ipSpecies long ipSpec = 0; for( int i = 0; i < nSpecies; i++) { ipSpec = i; if(dBaseStates[i][0].nelem()-1 == nelem && dBaseStates[i][0].IonStg()-1 == ion) break; ipSpec = -1; } if( ipSpec != -1) { for( int lvl = 0; lvl < dBaseSpecies[ipSpec].numLevels_max; lvl++) { //total level population totallevelpop += dBaseStates[ipSpec][lvl].Pop(); } } } for( int ion = 1; ion < nelem-1; ion++) { // Find ipSpecies long ipSpec = 0; for( int i = 0; i < nSpecies; i++) { ipSpec = i; if(dBaseStates[i][0].nelem()-1 == nelem && dBaseStates[i][0].IonStg()-1 == ion) break; ipSpec = -1; } //Energy Level Correction Factor - Scale energy levels to ground state of atomic Neon //using ionization potentials double ELCF = 0.; for( int i = 0; i < ion; i++ ) { ELCF += Heavy.Valence_IP_Ryd[nelem][i]*EVRYD; } if( ipSpec != -1) { //This requires a loop over energy levels for( int lvl = 0; lvl < dBaseSpecies[ipSpec].numLevels_max; lvl++) { //Give bound-free Destruction rates for the ground state, zero otherwise if( lvl == 0 && ion != nelem+1) { GcollBF = ionbal.CollIonRate_Ground[nelem][ion][0]; GphotoBF = PIR[ion]; GautoBF = autoion[ion]; } else { GcollBF = 0.; GphotoBF = 0.; GautoBF = 0.; } //Find the Destruction total rate Gtotal = dBaseStates[ipSpec][lvl].DestCollBB() + dBaseStates[ipSpec][lvl].DestPhotoBB() + GcollBF + GphotoBF + GautoBF; //Store Dest Rate Fractions double f_GcollBB, f_GphotoBB, f_GcollBF, f_GphotoBF, f_GautoBF; if( Gtotal != 0.) { f_GcollBB = dBaseStates[ipSpec][lvl].DestCollBB()/Gtotal; f_GphotoBB = dBaseStates[ipSpec][lvl].DestPhotoBB()/Gtotal; f_GcollBF = GcollBF/Gtotal; f_GphotoBF = GphotoBF/Gtotal; f_GautoBF = GautoBF/Gtotal; } else { f_GcollBB = 0.; f_GphotoBB = 0.; f_GcollBF = 0.; f_GphotoBF = 0.; f_GautoBF = 0.; } //Give bound-free Creation rates for the ground state, zero otherwise if( lvl == 0 && ion != 0 ) { QautoBF = dense.eden*ionbal.DR_Badnell_rate_coef[nelem][ion-1]; QphotoBF = dense.eden*ionbal.RR_rate_coef_used[nelem][ion-1]; QcollBF = dense.eden*ionbal.CotaRate[ion-1]; } else { QautoBF = 0.; QphotoBF = 0.; QcollBF = 0.; } //Find the Creation total rate Qtotal = dBaseStates[ipSpec][lvl].CreatCollBB() + 0.0 + QcollBF + QphotoBF + QautoBF; //Get Creation Rate Fractions double f_QcollBB, f_QphotoBB, f_QcollBF, f_QphotoBF, f_QautoBF; if( Qtotal != 0.) { f_QcollBB = dBaseStates[ipSpec][lvl].CreatCollBB()/Qtotal; f_QphotoBB = 0.0/Qtotal; f_QcollBF = QcollBF/Qtotal; f_QphotoBF = QphotoBF/Qtotal; f_QautoBF = QautoBF/Qtotal; } else { f_QcollBB = 0.; f_QphotoBB = 0.; f_QcollBF = 0.; f_QphotoBF = 0.; f_QautoBF = 0.; } //Section 4 print fprintf( save.ipPnunit[ipPun],"elev %li %3i %11.4e %11.6e %11.4e\n" " %11.3e %11.3e %11.3e %11.3e %11.3e %11.3e\n" " %11.3e %11.3e %11.3e %11.3e %11.3e %11.3e\n" " %i %i %i\n\n", nelem+1-ion, lvl+1, dBaseStates[ipSpec][lvl].g(), dBaseStates[ipSpec][lvl].energy().eV()+ELCF, dBaseStates[ipSpec][lvl].Pop()/totallevelpop, Gtotal*dBaseStates[ipSpec][lvl].Pop(), f_GcollBB, f_GphotoBB, f_GcollBF, f_GphotoBF, f_GautoBF, Qtotal, f_QcollBB, f_QphotoBB, f_QcollBF, f_QphotoBF, f_QautoBF, 0, 0, nouter[ion]); } } else printf("\nNot using :\t%li\t%i for levels\n",nelem+1,ion+1); } return; } # endif // implements the absorption option on the // set save line width command inline realnum PrettyTransmission(long j, realnum transmission) { if( save.ResolutionAbs < realnum(0.) ) // option to conserve energy return transmission; else { realnum corr = save.ResolutionAbs*rfield.widflx[j]/rfield.anu[j]; return realnum(max(0., 1. - (1.-transmission)*corr )); } } /*SaveResults1Line do single line of output for the save results and save line intensity commands */ /* the number of emission lines across one line of printout */ STATIC void SaveResults1Line( FILE * ioPUN, /* 4 char + null string */ const char *chLab, realnum wl, double xInten, const char *chFunction);/* null terminated string saying what to do */ /*SaveGaunts called by save gaunts command to output gaunt factors */ STATIC void SaveGaunts(FILE* ioPUN); /*SaveResults save results from save results command */ /*SaveResults1Line do single line of output for the save results and save line intensity commands */ STATIC void SaveResults(FILE* ioPUN); STATIC void SaveLineStuff( FILE * ioPUN, const char *chJob , realnum xLimit); /* save h emission, for chapter 4, routine is below */ STATIC void AGN_Hemis(FILE *ioPUN ); /*SaveNewContinuum produce the 'save new continuum' output */ STATIC void SaveNewContinuum(FILE * ioPUN ); /*SaveLineIntensity produce the 'save lines intensity' output */ STATIC void SaveLineIntensity(FILE * ioPUN , long int ipPun, realnum Threshold); char *chDummy; /*SaveFeII_cont return appropriate FeII banded continuum, inward, outward, or total */ STATIC realnum SaveFeII_cont( long int ipCont , long ipFeII_Cont_type ) { DEBUG_ENTRY( "SaveFeII_cont()" ); // [0]=wavelength, [1] inward; [2] outward, 3=> total if( ipFeII_Cont_type == 3 ) return( FeII_Cont[ipCont][1]+FeII_Cont[ipCont][2] ); else return(FeII_Cont[ipCont][ipFeII_Cont_type]); } void SaveDo( /* chTime is null terminated 4 char string, either "MIDL" or "LAST" */ const char *chTime) { bool lgTlkSav; long int ipPun, i, i1, ion, ipConMax, ipHi, ipLinMax, ipLo, ips, j, jj, limit, nelem; double ConMax, RateInter, av, conem, eps, flxatt, flxcor, flxin, flxref, flxtrn, fout, fref, fsum, opConSum, opLinSum, stage, sum, texc, xLinMax; DEBUG_ENTRY( "SaveDo()" ); /* * the "last" option on save command, to save on last iteration, * is parsed at the top of the loop in only one place. * no further action is needed at all for save last to work * ok throughout this routine */ /* * each branch can have a test whether chTime is or is not "LAST" * * if( strcmp(chTime,"LAST") == 0 ) <== print after iteration is complete * * if "LAST" then this is last call to routine after iteration complete * save only if "LAST" when results at end of iteration are needed * * if( strcmp(chTime,"LAST") != 0 ) <== print for every zone * * test for .not."LAST" is for every zone result, where you do not * want to save last zone twice */ /* return if no save to do */ if( save.nsave < 1 ) { return; } /* during a grid calculation this routine saves grid points after * cloudy is called. we may output it below */ if( grid.lgGrid ) { if( chTime[0]=='L' ) GridGatherInCloudy(); } // sort line labels if this is last call, this avoids multiple calls if several // output options need sorted labels and is safer since labels will be sorted in // case new code is added that reports the strong lines. The disadvantage is that // we sort even if the labels are not used if( strcmp(chTime,"LAST") == 0 ) { // sort emission line intensities so strongest lines are reported FindStrongestLineLabels(); } for( ipPun=0; ipPun < save.nsave; ipPun++ ) { if( save.lgPunHeader[ipPun] && !nMatch(save.chHeader[ipPun],save.chNONSENSE) ) { fprintf( save.ipPnunit[ipPun], "%s", save.chHeader[ipPun] ); save.lgPunHeader[ipPun] = false; } /* this global variable to remember where in the save stack we are */ save.ipConPun = ipPun; /* used to identify case where no key found */ bool lgNoHitFirstBranch = false; /* iterations.lgLastIt is true if this is last iteration * lgPunLstIter set true if 'last' key occurred on save command * normally is false. This will skip saving if last set and * this is not last iteration */ /* IMPORTANT: there is a second, identical if-statement halfway * down this routine. Any changes here should be copied there! */ if( iterations.lgLastIt || !save.lgPunLstIter[ipPun] || // if the sim aborted, make sure that the punch output is still done. // for MIDL output this is not very useful as all previous zones from // the iteration where the abort occured will be missing, but for LAST // output it is important to print at least placeholders in grid runs. ( lgAbort && strcmp(chTime,"LAST") == 0 ) || ( dynamics.lgTimeDependentStatic && dynamics.lgStatic_completed ) ) { if( strcmp(save.chSave[ipPun],"ABUN") == 0 ) { /* save abundances vs depth */ if( strcmp(chTime,"LAST") != 0 ) { fprintf( save.ipPnunit[ipPun], "%.2f", log10(MAX2(SMALLFLOAT,dense.gas_phase[ipHYDROGEN])) ); for( nelem=ipHELIUM; nelem < LIMELM; nelem++ ) { /* >>chng 05 feb 03, protect against non-positive abundances, * bug caught by Marcelo Castellanos */ fprintf( save.ipPnunit[ipPun], "\t%.2f", log10(MAX2(SMALLFLOAT,dense.gas_phase[nelem])) ); } fprintf( save.ipPnunit[ipPun], "\n" ); } } else if( strcmp(save.chSave[ipPun],"21CM") == 0 ) { /* save information about 21 cm line */ if( strcmp(chTime,"LAST") != 0 ) { fprintf( save.ipPnunit[ipPun], "%.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", /* depth, cm */ radius.depth_mid_zone, hyperfine.Tspin21cm , phycon.te , /* temperature from Lya - 21 cm optical depth ratio */ 3.84e-7* iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauCon() / SDIV( HFLines[0].Emis().TauCon() ), /*TexcLine( &iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) ),*/ (*HFLines[0].Lo()).Pop() , (*HFLines[0].Hi()).Pop() , OccupationNumberLine( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) ), HFLines[0].Emis().TauCon() , iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauCon(), HFLines[0].Emis().PopOpc(), /* term in () is density (cm-3) of 1s, so this is n(1s) / Ts */ (iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop())/SDIV( hyperfine.Tspin21cm), /* why was above multiplied by this following term? */ /* *HFLines[0].EnergyErg/BOLTZMANN/4.,*/ HFLines[0].Emis().TauIn(), TexcLine( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) ) , colden.H0_ov_Tspin, /*>>chng 27 mar, GS, integrated 21cm spin temperature*/ colden.H0_21cm_lower, colden.H0_21cm_upper, -0.068/log((colden.H0_21cm_upper/3.)/colden.H0_21cm_lower) ); } } else if( strcmp(save.chSave[ipPun],"AGES") == 0 ) { /* save timescales vs depth */ if( strcmp(chTime,"LAST") != 0 ) { int ipCO, ipOH; ipCO = findspecies("CO")->index; ipOH = findspecies("OH")->index; fprintf( save.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", /* depth, cm */ radius.depth_mid_zone, /* cooling timescale */ dense.pden*BOLTZMANN*1.5*phycon.te/ thermal.htot, /* H2 destruction timescale */ timesc.time_H2_Dest_here, /* CO destruction timescale */ 1./SDIV((ipCO != -1) ? mole.species[ipCO].snk : 0.), /* OH destruction timescale */ 1./SDIV((ipOH != -1) ? mole.species[ipOH].snk : 0.), /* H recombination timescale */ 1./(dense.eden*2.90e-10/(phycon.te70*phycon.te10/phycon.te03)) ); } } else if( strcmp(save.chSave[ipPun]," AGN") == 0 ) { if( strcmp(chTime,"LAST") == 0 ) { if( strcmp( save.chSaveArgs[ipPun], "HECS" ) == 0 ) { /* this routine is in helike.c */ AGN_He1_CS(save.ipPnunit[ipPun]); } if( strcmp( save.chSaveArgs[ipPun], "HEMI" ) == 0 ) { /* save h emiss, for chapter 4, routine is below */ AGN_Hemis(save.ipPnunit[ipPun]); } else { fprintf( ioQQQ, " SaveDo does not recognize flag %4.4s set for AGN save. This is impossible.\n", save.chSave[ipPun] ); ShowMe(); cdEXIT(EXIT_FAILURE); } } } else if( strcmp(save.chSave[ipPun],"MONI") == 0 ) { if( strcmp(chTime,"LAST") == 0 ) { /* save the monitor output */ lgCheckMonitors(save.ipPnunit[ipPun]); } } else if( strcmp(save.chSave[ipPun],"AVER") == 0 ) { if( strcmp(chTime,"LAST") == 0 ) { /* save the averages output */ save_average( ipPun ); } } else if( strncmp(save.chSave[ipPun],"CHA",3) == 0 ) { if( strcmp(chTime,"LAST") == 0 ) { /* one of the charge transfer options, all in chargtran.c */ ChargTranPun( save.ipPnunit[ipPun] , save.chSave[ipPun] ); } } else if( strcmp( save.chSave[ipPun],"CHIA") == 0) { static bool lgRunOnce = true; if( lgRunOnce ) { lgRunOnce = false; // save chianti collision data in physical units int ipLo = 0; int ipHi = 0; double fupsilon = 0.; double initTemp = 4.0; double finalTemp = 9.1; double stepTemp = 0.2; fprintf( save.ipPnunit[ipPun],"Species\tLo\tHi\tWLAng"); for(double logtemp = initTemp;logtemp < finalTemp;logtemp = logtemp + stepTemp ) { fprintf( save.ipPnunit[ipPun],"\t%2.1f",logtemp); } fprintf( save.ipPnunit[ipPun],"\n"); for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies) { for( EmissionList::iterator tr=dBaseTrans[ipSpecies].Emis().begin(); tr != dBaseTrans[ipSpecies].Emis().end(); ++tr) { if( dBaseTrans[ipSpecies].chLabel() == "Chianti" ) { ipLo = tr->Tran().ipLo(); ipHi = tr->Tran().ipHi(); fprintf( save.ipPnunit[ipPun],"%s\t%i\t%i\t", dBaseSpecies[ipSpecies].chLabel,ipLo+1,ipHi+1); fprintf( save.ipPnunit[ipPun],"%5.3e",tr->Tran().WLAng()); for(double logtemp = initTemp;logtemp < finalTemp;logtemp = logtemp + stepTemp ) { fupsilon = CHIANTI_Upsilon(ipSpecies, ipELECTRON, ipHi, ipLo, pow(10.,logtemp)); fprintf( save.ipPnunit[ipPun],"\t%5.3e",fupsilon); } fprintf( save.ipPnunit[ipPun],"\n"); } } } } } else if( strcmp(save.chSave[ipPun],"COOL") == 0 || strcmp(save.chSave[ipPun],"EACH") == 0 ) { /* save cooling, routine in file of same name */ if( strcmp(chTime,"LAST") != 0 ) CoolSave(save.ipPnunit[ipPun], save.chSave[ipPun]); } else if( strcmp(save.chSave[ipPun],"DOMI") == 0 ) { /* save dominant rates */ if( strcmp(chTime,"LAST") != 0 ) { molecule *debug_species = findspecies( save.chSpeciesDominantRates[ipPun] ); if (debug_species == null_mole) { fprintf( ioQQQ,"Error in SAVE DOMINANT RATES, species %s not found\n", save.chSpeciesDominantRates[ipPun]); } else { fprintf( save.ipPnunit[ipPun], "%e\t%e\t", radius.depth_mid_zone, mole.species[ debug_species->index ].column ); mole_dominant_rates( debug_species, save.ipPnunit[ipPun] ); } } } else if( strncmp(save.chSave[ipPun],"DYN" , 3) == 0 ) { /* save dynamics xxx, information about dynamical solutions */ if( strcmp(chTime,"LAST") != 0 ) DynaSave(save.ipPnunit[ipPun] ,save.chSave[ipPun][3] ); } else if( strcmp(save.chSave[ipPun],"ENTH") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) fprintf( save.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", radius.depth_mid_zone, phycon.EnthalpyDensity, phycon.EnergyExcitation, phycon.EnergyIonization, phycon.EnergyBinding , 0.5*POW2(wind.windv)*dense.xMassDensity , /* KE */ 5./2.*pressure.PresGasCurr , /* thermal plus PdV work */ magnetic.EnthalpyDensity); /* magnetic terms */ } else if( strcmp(save.chSave[ipPun],"COLU") == 0 ) { /* save column densities */ if( strcmp(chTime,"LAST") == 0 ) { PrtColumns(save.ipPnunit[ipPun],"TABLE",ipPun); } } else if( strcmp(save.chSave[ipPun],"COLS") == 0 ) { /* save some column densities */ if( strcmp(chTime,"LAST") == 0 ) { /*TODO unify following with PrtColumns above */ save_colden( save.ipPnunit[ipPun] ); } } else if( strcmp(save.chSave[ipPun],"COMP") == 0 ) { /* Compton energy exchange coefficients */ if( strcmp(chTime,"LAST") != 0 ) { for( jj=0; jj>chng 06 apr 03, add every option to do every zone */ /* if save every is set then nSaveEveryZone must be positive * was init to -1 */ bool lgPrintThis =false; if( save.lgSaveEveryZone[ipPun] ) { /* this branch, every option is on line so want to print every n zone */ if( strcmp(chTime,"LAST") != 0 ) { /* not last zone - print first and intermediate cases */ if( nzone==1 ) { lgPrintThis = true; } else if( nzone%save.nSaveEveryZone[ipPun]==0 ) { lgPrintThis = true; } } else { /* this is last zone, print only if did not trip on above */ if( nzone!=1 && nzone%save.nSaveEveryZone[ipPun]!=0 ) { lgPrintThis = true; } } } else { /* this branch, not "every", so only print the last zone */ if( strcmp(chTime,"LAST") == 0 ) lgPrintThis = true; } ASSERT( !save.lgSaveEveryZone[ipPun] || save.nSaveEveryZone[ipPun]>0 ); if( lgPrintThis ) { if( save.lgSaveEveryZone[ipPun] && nzone!=1) fprintf( save.ipPnunit[ipPun], "%s\n", save.chHashString ); /* option to also print same arrays but for cumulative arrays */ int nEmType = (int)save.punarg[ipPun][0]; ASSERT( nEmType==0 || nEmType==1 ); const realnum *trans_coef_total=rfield.getCoarseTransCoef(); for( j=0; j>chng 96 oct 22, format of anu to .5 to resolve energy mesh near 1 Ryd */ fprintf( save.ipPnunit[ipPun], "%.5e\t%.3e\n", AnuUnit(rfield.AnuOrg[j]), flxin ); } } } else if( strcmp(save.chSave[ipPun],"CONG") == 0 ) { /* save emitted grain continuum in optically thin limit */ if( strcmp(chTime,"LAST") == 0 ) { /* GrainMakeDiffuse broke out emission into types * according to matType */ for( j=0; j>>>>>>>>>>>>\n Reflected continuum not predicted when SPHERE is set.\n" ); cdEXIT(EXIT_FAILURE); } for( j=0; j 1e-25 ) { av = rfield.ConRefIncid[0][j]/rfield.flux_total_incident[0][j]; } else { av = 0.; } /* >>chng 96 oct 22, format of anu to .5 to resolve energy mesh near 1 Ryd */ fprintf( save.ipPnunit[ipPun], "%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4s\n", AnuUnit(rfield.AnuOrg[j]), flxref, fref, flxref + fref, av, rfield.chContLabel[j] ); } } } else if( strcmp(save.chSave[ipPun],"CNVE") == 0 ) { /* the save convergence error command */ if( strcmp(chTime,"LAST") != 0 ) { fprintf( save.ipPnunit[ipPun], "%.5e\t%li\t%.4e\t%.4f\t%.4e\t%.4e\t%.3f\t%.4e\t%.4e\t%.4f\n", radius.depth_mid_zone, conv.nPres2Ioniz, pressure.PresTotlCurr, pressure.PresTotlError*100., dense.EdenTrue, dense.eden, (dense.EdenTrue - dense.eden)*100./dense.EdenTrue, thermal.htot, thermal.ctot, (thermal.htot - thermal.ctot)*100./thermal.htot ); } } else if( strcmp(save.chSave[ipPun],"CONB") == 0 ) { /* save continuum bins binning */ /* set pointer for possible change in units of energy in continuum * AnuUnit will give anu in whatever units were set with save units */ if( strcmp(chTime,"LAST") != 0 ) { for( j=0; j>chng 96 oct 22, format of anu to .5e to resolve energy mesh near 1 Ryd */ fprintf( save.ipPnunit[ipPun], "%.5e\t%.3e\t%.3e\t%.3e\t%4.4s\t%4.4s\n", AnuUnit(rfield.AnuOrg[j]), flxref, conem, flxref + conem, rfield.chLineLabel[j], rfield.chContLabel[j] ); } } } /* save fine continuum command */ else if( strcmp(save.chSave[ipPun],"CONf") == 0 ) { if( strcmp(chTime,"LAST") == 0 ) { long nu_hi , nskip; if( save.punarg[ipPun][0] > 0. ) /* possible lower bounds to energy range - * 0 if not set with range option*/ j = ipFineCont( save.punarg[ipPun][0] ); else j = 0; /* upper limit set with range option */ if( save.punarg[ipPun][1]> 0. ) nu_hi = ipFineCont( save.punarg[ipPun][1]); else nu_hi = rfield.nfine; /* number of cells to bring together, default is 10 */ nskip = (long)save.punarg[ipPun][2]; do { realnum sum1 = rfield.fine_opt_depth[j]; realnum xnu = rfield.fine_anu[j]; for( jj=1; jj 0. ) sum = 1./sum; else sum = 1.; fsum = 0.; for( j=i1-1; j>chng 01 nov 22, format to c-friendly */ if( (RateInter >= save.punarg[ipPun][1]) && (flxcor > SMALLFLOAT) ) { lgPrt = true; /* >>chng 96 oct 22, format of anu to 11.5 to resolve energy mesh near 1 Ryd */ fprintf( save.ipPnunit[ipPun], "%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", j, AnuUnit(rfield.AnuOrg[j]), flxcor, flxcor*opac.opacity_abs[j], rfield.flux[0][j]/flxcor, rfield.otslin[j]/flxcor, rfield.otscon[j]/flxcor, (rfield.outlin[0][j] + rfield.outlin_noplot[j])/flxcor, rfield.ConInterOut[j]/flxcor, RateInter, fsum*sum, rfield.chLineLabel[j], rfield.chContLabel[j] ); } } if( !lgPrt ) { /* entered logical block but did not print anything */ fprintf(save.ipPnunit[ipPun], " punchdo, the PUNCH IONIZING CONTINUUM command " "did not find a strong point, sum and fsum were %.2e %.2e\n", sum,fsum); fprintf(save.ipPnunit[ipPun], " punchdo, the low-frequency energy was %.5e Ryd\n", rfield.anu[i1-1]); fprintf(save.ipPnunit[ipPun], " You can reset the threshold for the lowest fractional " "interaction to print with the second number of the save command\n" " The fraction was %.3f and this was too large.\n", save.punarg[ipPun][1]); } } } #ifdef USE_NLTE7 else if( strcmp(save.chSave[ipPun],"CONl") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) { int nEmType = (int)save.punarg[ipPun][0]; ASSERT( nEmType==0 || nEmType==1 ); for( j=0; j= lowlim && NRGeV <= highlim ) { /* photon energy in appropriate energy or wavelength units */ fprintf( save.ipPnunit[ipPun],"%14.8e ", NRGeV ); /* print zeroes for BB BF and FF */ fprintf( save.ipPnunit[ipPun],"%14.8e %14.8e %14.8e ",0.,0.,0.); /* total cont */ fprintf( save.ipPnunit[ipPun],"%14.8e\n", (flxref + flxtrn)/NRGeV ); } } } } #endif else if( strcmp(save.chSave[ipPun],"CONS") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) { // continuum volume emissivity and opacity as a function of radius // command was "save continuum emissivity" */ if( save.ipEmisFreq[ipPun] < 0 ) save.ipEmisFreq[ipPun] = ipoint(save.emisfreq[ipPun].Ryd()); j = save.ipEmisFreq[ipPun]-1; fprintf( save.ipPnunit[ipPun], "%.14e\t%.14e\t%.5e\t%.5e\t%.5e\n", radius.Radius_mid_zone, radius.depth_mid_zone, rfield.anu2[j]*rfield.ConEmitLocal[nzone][j]/rfield.widflx[j]*EN1RYD, opac.opacity_abs[j], opac.opacity_sct[j] ); } } else if( strcmp(save.chSave[ipPun],"CORA") == 0 ) { /* save raw continuum */ /* set pointer for possible change in units of energy in continuum * AnuUnit will give anu in whatever units were set with save units */ if( strcmp(chTime,"LAST") == 0 ) { /* this option to save all raw ionizing continuum */ for( j=0;j>chng 97 jul 10, remove SaveLWidth from this one only since * we must conserve energy even in lines * >>chng 07 apr 26 include transmission coefficient */ flxatt = rfield.flux[0][j]*rfield.anu2[j]*EN1RYD/ rfield.widflx[j]*radius.r1r0sq*trans_coef_total[j]; /*conem = (rfield.ConOutNoInter[j] + rfield.ConInterOut[j]+rfield.outlin[0][j])* rfield.anu2[j]; conem *= radius.r1r0sq*EN1RYD*geometry.covgeo;*/ /* >>chng 00 jan 03, above did not include all contributors. * Pasted in below from usual * save continuum command */ /* >>chng 04 jul 15, removed factor of save.SaveLWidth - * this should not be there to conserve energy, as explained in hazy * where command was documented, and in comment above. caught by PvH */ /* >>chng 04 jul 23, incorrect use of outlin - before multiplied by an2, * quantity should be photons per Ryd, since init quantity is * photons per cell. Must div by widflx. caught by PvH */ /*conem = (rfield.ConEmitOut[0][j]/rfield.widflx[j]*rfield.anu2[j] + rfield.outlin[0][j]*rfield.anu[j])*radius.r1r0sq* EN1RYD*geometry.covgeo;*/ conem = (rfield.ConEmitOut[0][j] + rfield.outlin[0][j]) / rfield.widflx[j]* rfield.anu2[j]*radius.r1r0sq*EN1RYD*geometry.covgeo; flxtrn = conem + flxatt; /* use AnuOrg here instead of anu since probably * going to be used by table read * and we want good anu array for sanity check*/ fprintf( save.ipPnunit[ipPun],"%.5e\t%.3e\t%.3e\n", AnuUnit(rfield.AnuOrg[j]), flxtrn, trans_coef_total[j] ); } } } else if( strcmp(save.chSave[ipPun],"CON2") == 0 ) { /* save total two-pohoton continuum */ if( strcmp(chTime,"LAST") == 0 ) { /* this option to save diffuse continuum */ for( j=0; jpure_sc1[j]*gv.bin[nd]->dstAbund; } /* finally, scattering including effects of forward scattering */ fprintf( save.ipPnunit[ipPun], "%.2e\t", scat ); fprintf( save.ipPnunit[ipPun], "%.2e\n", gv.dstsc[j]/SDIV(gv.dstab[j] + gv.dstsc[j]) ); } } } /* save grain abundance and save grain D/G ratio commands */ else if( strcmp(save.chSave[ipPun],"DUSA") == 0 || strcmp(save.chSave[ipPun],"DUSD") == 0 ) { bool lgDGRatio = ( strcmp(save.chSave[ipPun],"DUSD") == 0 ); /* grain abundance */ if( strcmp(chTime,"LAST") != 0 ) { /* used to print header exactly one time */ static bool lgMustPrtHeaderDRRatio = true, lgMustPrtHeaderAbundance=true; /* print grain header first if this has not yet been done */ if( ( lgMustPrtHeaderDRRatio && lgDGRatio ) || ( lgMustPrtHeaderAbundance && !lgDGRatio ) ) { /* only print one header for each case, but must * track separately if both used in same sim */ if( lgMustPrtHeaderDRRatio && lgDGRatio ) lgMustPrtHeaderDRRatio = false; else if( lgMustPrtHeaderAbundance &&!lgDGRatio ) lgMustPrtHeaderAbundance = false; fprintf( save.ipPnunit[ipPun], "#Depth" ); for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab ); fprintf( save.ipPnunit[ipPun], "\ttotal\n" ); /* now print grain radius converting from cm to microns */ fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" ); for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 ); fprintf( save.ipPnunit[ipPun], "\txxxx\n" ); } fprintf( save.ipPnunit[ipPun], " %.5e", radius.depth_mid_zone ); /* grain abundance per bin in g/cm^3 */ double total = 0.; for( size_t nd=0; nd < gv.bin.size(); ++nd ) { double abund = gv.bin[nd]->IntVol*gv.bin[nd]->dustp[0]* gv.bin[nd]->cnv_H_pCM3; if( lgDGRatio ) abund /= dense.xMassDensity; fprintf( save.ipPnunit[ipPun], "\t%.3e", abund ); total += abund; } fprintf( save.ipPnunit[ipPun], "\t%.3e\n", total ); } } else if( strcmp(save.chSave[ipPun],"DUSP") == 0 ) { /* grain potential */ if( strcmp(chTime,"LAST") != 0 ) { /* used to print header exactly one time */ static bool lgMustPrtHeader = true; /* do labels first if this is first zone */ if( lgMustPrtHeader ) { /* first print string giving grain id */ fprintf( save.ipPnunit[ipPun], "#Depth" ); for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab ); fprintf( save.ipPnunit[ipPun], "\n" ); /* now print grain radius converting from cm to microns */ fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" ); for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 ); fprintf( save.ipPnunit[ipPun], "\n" ); /* don't need to do this, ever again */ lgMustPrtHeader = false; } fprintf( save.ipPnunit[ipPun], " %.5e", radius.depth_mid_zone ); /* grain potential in eV */ for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->dstpot*EVRYD ); fprintf( save.ipPnunit[ipPun], "\n" ); } } else if( strcmp(save.chSave[ipPun],"DUSR") == 0 ) { /* grain H2 formation rates */ if( strcmp(chTime,"LAST") != 0 ) { /* used to print header exactly one time */ static bool lgMustPrtHeader = true; /* do labels first if this is first zone */ if( lgMustPrtHeader ) { /* first print string giving grain id */ fprintf( save.ipPnunit[ipPun], "#Depth" ); for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab ); fprintf( save.ipPnunit[ipPun], "\n" ); /* now print grain radius converting from cm to microns */ fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" ); for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 ); fprintf( save.ipPnunit[ipPun], "\n" ); /* don't need to do this, ever again */ lgMustPrtHeader = false; } fprintf( save.ipPnunit[ipPun], " %.5e", radius.depth_mid_zone ); /* grain formation rate for H2 */ for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->rate_h2_form_grains_used ); fprintf( save.ipPnunit[ipPun], "\n" ); } } else if( strcmp(save.chSave[ipPun],"DUST") == 0 ) { /* grain temperatures - K*/ if( strcmp(chTime,"LAST") != 0 ) { /* used to print header exactly one time */ static bool lgMustPrtHeader = true; /* do labels first if this is first zone */ if( lgMustPrtHeader ) { /* first print string giving grain id */ fprintf( save.ipPnunit[ipPun], "#Depth" ); for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab ); fprintf( save.ipPnunit[ipPun], "\n" ); /* now print grain radius converting from cm to microns */ fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" ); for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 ); fprintf( save.ipPnunit[ipPun], "\n" ); /* don't need to do this, ever again */ lgMustPrtHeader = false; } fprintf( save.ipPnunit[ipPun], " %.5e", radius.depth_mid_zone ); for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->tedust ); fprintf( save.ipPnunit[ipPun], "\n" ); } } else if( strcmp(save.chSave[ipPun],"DUSC") == 0 ) { /* save grain charge - eden from grains and * charge per grain in electrons / grain */ if( strcmp(chTime,"LAST") != 0 ) { /* used to print header exactly one time */ static bool lgMustPrtHeader = true; /* do labels first if this is first zone */ if( lgMustPrtHeader ) { /* first print string giving grain id */ fprintf( save.ipPnunit[ipPun], "#Depth\tne(grn)" ); for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab ); fprintf( save.ipPnunit[ipPun], "\n" ); /* now print grain radius converting from cm to microns */ fprintf( save.ipPnunit[ipPun], "#grn rad (mic)\tne\t" ); for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 ); fprintf( save.ipPnunit[ipPun], "\n" ); /* don't need to do this, ever again */ lgMustPrtHeader = false; } fprintf( save.ipPnunit[ipPun], " %.5e\t%.4e", radius.depth_mid_zone , /* electron density contributed by grains, in e/cm^3, * positive number means grain supplied free electrons */ gv.TotalEden ); /* average charge per grain in electrons */ for( size_t nd=0; nd < gv.bin.size(); ++nd ) { fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AveDustZ ); } fprintf( save.ipPnunit[ipPun], "\n" ); } } else if( strcmp(save.chSave[ipPun],"DUSH") == 0 ) { /* grain heating */ if( strcmp(chTime,"LAST") != 0 ) { /* used to print header exactly one time */ static bool lgMustPrtHeader = true; /* save grain charge, but do labels first if this is first zone */ if( lgMustPrtHeader ) { /* first print string giving grain id */ fprintf( save.ipPnunit[ipPun], "#Depth" ); for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab ); fprintf( save.ipPnunit[ipPun], "\n" ); /* now print grain radius converting from cm to microns */ fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" ); for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 ); fprintf( save.ipPnunit[ipPun], "\n" ); /* don't need to do this, ever again */ lgMustPrtHeader = false; } fprintf( save.ipPnunit[ipPun], " %.5e", radius.depth_mid_zone ); /* grain heating */ for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->GasHeatPhotoEl ); fprintf( save.ipPnunit[ipPun], "\n" ); } } else if( strcmp(save.chSave[ipPun],"DUSV") == 0 ) { /* grain drift velocities */ if( strcmp(chTime,"LAST") != 0 ) { /* used to print header exactly one time */ static bool lgMustPrtHeader = true; /* save grain velocity, but do labels first if this is first zone */ if( lgMustPrtHeader ) { /* first print string giving grain id */ fprintf( save.ipPnunit[ipPun], "#Depth" ); for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab ); fprintf( save.ipPnunit[ipPun], "\n" ); /* now print grain radius converting from cm to microns */ fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" ); for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 ); fprintf( save.ipPnunit[ipPun], "\n" ); /* don't need to do this, ever again */ lgMustPrtHeader = false; } fprintf( save.ipPnunit[ipPun], " %.5e", radius.depth_mid_zone ); /* grain drift velocity in km/s */ for( size_t nd=0; nd < gv.bin.size(); ++nd ) fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->DustDftVel*1e-5 ); fprintf( save.ipPnunit[ipPun], "\n" ); } } /* >>chng 02 dec 30, separated scattering cross section and asymmetry factor, PvH */ else if( strcmp(save.chSave[ipPun],"DUSQ") == 0 ) { /* save grain Qs */ if( strcmp(chTime,"LAST") == 0 ) { for( j=0; j < rfield.nflux; j++ ) { fprintf( save.ipPnunit[ipPun], "%11.4e", rfield.anu[j] ); for( size_t nd=0; nd < gv.bin.size(); nd++ ) { fprintf( save.ipPnunit[ipPun], "%9.1e%9.1e", gv.bin[nd]->dstab1[j]*4./gv.bin[nd]->IntArea, gv.bin[nd]->pure_sc1[j]*gv.bin[nd]->asym[j]*4./gv.bin[nd]->IntArea ); } fprintf( save.ipPnunit[ipPun], "\n" ); } } } else if( strcmp(save.chSave[ipPun],"ELEM") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) { realnum renorm = 1.f; /* this is the index for the atomic number on the physical scale */ /* >>chng 04 nov 23, use c scale throughout */ nelem = (long int)save.punarg[ipPun][0]; ASSERT( nelem >= ipHYDROGEN ); /* don't do this if element is not turned on */ if( dense.lgElmtOn[nelem] ) { /* >>chng 04 nov 23, add density option, leave as cm-3 * default is still norm to total of that element */ if( save.punarg[ipPun][1] == 0 ) renorm = dense.gas_phase[nelem]; fprintf( save.ipPnunit[ipPun], " %.5e", radius.depth_mid_zone ); for( j=0; j <= (nelem + 1); ++j) { fprintf( save.ipPnunit[ipPun], "\t%.2e", dense.xIonDense[nelem][j]/renorm ); } if( nelem==ipHYDROGEN ) { /* H2 */ fprintf( save.ipPnunit[ipPun], "\t%.2e", hmi.H2_total/renorm ); } /* >>chng 04 nov 23 add C and O fine structure pops */ else if( nelem==ipCARBON ) { fprintf( save.ipPnunit[ipPun], "\t%.2e\t%.2e\t%.2e", colden.C1Pops[0]/renorm, colden.C1Pops[1]/renorm, colden.C1Pops[2]/renorm); fprintf( save.ipPnunit[ipPun], "\t%.2e\t%.2e", colden.C2Pops[0]/renorm, colden.C2Pops[1]/renorm); fprintf( save.ipPnunit[ipPun], "\t%.2e", findspecieslocal("CO")->den/renorm ); } else if( nelem==ipOXYGEN ) { fprintf( save.ipPnunit[ipPun], "\t%.2e\t%.2e\t%.2e", colden.O1Pops[0]/renorm, colden.O1Pops[1]/renorm, colden.O1Pops[2]/renorm); } fprintf( save.ipPnunit[ipPun], "\n" ); } } } else if( strcmp(save.chSave[ipPun],"RECA") == 0 ) { /* this will create table for AGN3 then exit, * routine is in makerecom.c */ ion_recombAGN( save.ipPnunit[ipPun] ); cdEXIT(EXIT_FAILURE); } else if( strcmp(save.chSave[ipPun],"RECE") == 0 ) { /* save recombination efficiencies, * option turned on with the "save recombination efficiencies" command * output for the save recombination coefficients command is actually * produced by a series of routines, as they generate the recombination * coefficients. these include * dielsupres, helium, hydrorecom, iibod, and makerecomb*/ fprintf( save.ipPnunit[ipPun], "%12.4e %12.4e %12.4e %12.4e\n", iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].RadRecomb[ipRecRad], iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].RadRecomb[ipRecNetEsc] , iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].RadRecomb[ipRecRad], iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].RadRecomb[ipRecNetEsc]); } else if( strncmp(save.chSave[ipPun],"FEc" , 3 ) == 0 ) { /* FeII continuum */ if( strcmp(chTime,"LAST") == 0 ) { if( save.punarg[ipPun][0] == 1 ) { // row format used for grids - one time print of wavelength grid first char chTempHeader[100]; long ipFeII_Cont_type; if( strcmp(save.chSave[ipPun],"FEcI") == 0 ) { strcpy(chTempHeader , "#FeII inward: Wl(A)\tInt[erg cm-2 s-1]\n"); ipFeII_Cont_type = 1; } else if( strcmp(save.chSave[ipPun],"FEcO") == 0 ) { strcpy(chTempHeader , "#FeII outward: Wl(A)\tInt[erg cm-2 s-1]\n"); ipFeII_Cont_type = 2; } else if( strcmp(save.chSave[ipPun],"FEcT") == 0 ) { strcpy(chTempHeader , "#FeII total: Wl(A)\tInt[erg cm-2 s-1]\n"); ipFeII_Cont_type = 3; } else TotalInsanity(); if( save.lgPunHeader[ipPun] && (!grid.lgGrid || optimize.nOptimiz==0) ) { // one time print of FeII continuum header j = 0; fprintf( save.ipPnunit[ipPun], "%s%.5f", chTempHeader , AnuUnit((realnum)RYDLAM/((FeII_Cont[j][0]+FeII_Cont[j+1][0])/2.f) )); for( j=1; j < nFeIIConBins; j++ ) { fprintf( save.ipPnunit[ipPun], "\t%.5e", AnuUnit((realnum)RYDLAM/((FeII_Cont[j][0]+FeII_Cont[j+1][0])/2.f) )); } fprintf( save.ipPnunit[ipPun], "\n"); save.lgPunHeader[ipPun] = false; } // give intensities fprintf( save.ipPnunit[ipPun], "%.2f", SaveFeII_cont( 0 , ipFeII_Cont_type )); for( j=1; j < nFeIIConBins; j++ ) { fprintf( save.ipPnunit[ipPun], "\t%e", SaveFeII_cont( j , ipFeII_Cont_type ) ); } fprintf( save.ipPnunit[ipPun], "\n"); } else if( save.punarg[ipPun][0] == 2 ) { // the default, four columns, wl, total, inward, outward // one time print of header if( save.lgPunHeader[ipPun] && (!grid.lgGrid || optimize.nOptimiz==0) ) { fprintf( save.ipPnunit[ipPun], "FeiI wl, total, inward, outward\n"); save.lgPunHeader[ipPun] = false; } for( j=0; j < nFeIIConBins; j++ ) { fprintf( save.ipPnunit[ipPun], "%.5f", AnuUnit((realnum)RYDLAM/((FeII_Cont[j][0]+FeII_Cont[j+1][0])/2.f) )); fprintf( save.ipPnunit[ipPun], "\t%e", SaveFeII_cont( j , 3 ) );// total fprintf( save.ipPnunit[ipPun], "\t%e", SaveFeII_cont( j , 1 ) );// inward fprintf( save.ipPnunit[ipPun], "\t%e", SaveFeII_cont( j , 2 ) );// outward fprintf( save.ipPnunit[ipPun], "\n"); } } } } /* save column densities */ else if( strcmp(save.chSave[ipPun],"FENl") == 0 ) { if( strcmp(chTime,"LAST") == 0 ) { /* save FeII level energies and stat weights, followed by column density */ FeIIPunchColden( save.ipPnunit[ipPun] ); } } else if( strcmp(save.chSave[ipPun],"FE2l") == 0 ) { if( strcmp(chTime,"LAST") == 0 ) { /* save FeII level energies and stat weights */ FeIIPunchLevels( save.ipPnunit[ipPun] ); } } else if( strcmp(save.chSave[ipPun],"FEli") == 0 ) { if( strcmp(chTime,"LAST") == 0 ) { /* save line intensities, routine is in atom_feii.c */ FeIISaveLines( save.ipPnunit[ipPun] ); } } else if( strcmp(save.chSave[ipPun],"FE2o") == 0 ) { if( strcmp(chTime,"LAST") == 0 ) { /* save line optical depths, routine is in atom_feii.c */ FeIIPunchOpticalDepth( save.ipPnunit[ipPun] ); } } else if( strcmp(save.chSave[ipPun],"FRED") == 0 ) { /* set with save Fred command, this punches some stuff from * Fred Hamann's dynamics project */ if( strcmp(chTime,"LAST") != 0 ) { /* Fred's list */ fprintf( save.ipPnunit[ipPun], "%.5e\t%.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\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\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", radius.Radius, radius.depth ,wind.windv/1e5, wind.dvdr, dense.gas_phase[ipHYDROGEN], dense.eden , phycon.te, wind.AccelLine , wind.AccelCont , wind.fmul , // acceleration in this zone due to electron scattering, // if incident SED was not attenuated pressure.pinzon_PresIntegElecThin/dense.xMassDensity/radius.drad_x_fillfac , mean.xIonMean[0][ipHYDROGEN][0][0] , mean.xIonMean[0][ipHYDROGEN][1][0] , mean.xIonMean[0][ipHELIUM][0][0] , mean.xIonMean[0][ipHELIUM][1][0] , mean.xIonMean[0][ipHELIUM][2][0] , mean.xIonMean[0][ipCARBON][1][0] , mean.xIonMean[0][ipCARBON][2][0] , mean.xIonMean[0][ipCARBON][3][0] , mean.xIonMean[0][ipOXYGEN][0][0] , mean.xIonMean[0][ipOXYGEN][1][0] , mean.xIonMean[0][ipOXYGEN][2][0] , mean.xIonMean[0][ipOXYGEN][3][0] , mean.xIonMean[0][ipOXYGEN][4][0] , mean.xIonMean[0][ipOXYGEN][5][0] , mean.xIonMean[0][ipOXYGEN][6][0] , mean.xIonMean[0][ipOXYGEN][7][0] , dense.xIonDense[ipHYDROGEN][0] , dense.xIonDense[ipHYDROGEN][1] , dense.xIonDense[ipHELIUM][0] , dense.xIonDense[ipHELIUM][1] , dense.xIonDense[ipHELIUM][2] , dense.xIonDense[ipCARBON][1] , dense.xIonDense[ipCARBON][2] , dense.xIonDense[ipCARBON][3] , dense.xIonDense[ipOXYGEN][0] , dense.xIonDense[ipOXYGEN][1] , dense.xIonDense[ipOXYGEN][2] , dense.xIonDense[ipOXYGEN][3] , dense.xIonDense[ipOXYGEN][4] , dense.xIonDense[ipOXYGEN][5] , dense.xIonDense[ipOXYGEN][6] , dense.xIonDense[ipOXYGEN][7] , mean.xIonMean[0][ipMAGNESIUM][1][0] , dense.xIonDense[ipMAGNESIUM][1], TauLines[ipT1032].Emis().TauIn() , TauLines[ipT1032].Emis().TauCon() ); } } else if( strcmp(save.chSave[ipPun],"FE2d") == 0 ) { /* save some departure coefficients for large FeII atom */ if( strcmp(chTime,"LAST") != 0 ) FeIIPunDepart(save.ipPnunit[ipPun],false); } else if( strcmp(save.chSave[ipPun],"FE2D") == 0 ) { /* save all departure coefficients for large FeII atom */ if( strcmp(chTime,"LAST") != 0 ) FeIIPunDepart(save.ipPnunit[ipPun],true); } else if( strcmp(save.chSave[ipPun],"FE2p") == 0 ) { bool lgFlag = false; if( save.punarg[ipPun][2] ) lgFlag = true; /* save small subset of level populations for large FeII atom */ if( strcmp(chTime,"LAST") != 0 ) FeIIPunPop(save.ipPnunit[ipPun],false,0,0, lgFlag); } else if( strcmp(save.chSave[ipPun],"FE2P") == 0 ) { bool lgFlag = false; if( save.punarg[ipPun][2] ) lgFlag = true; /* save range of level populations for large FeII atom */ if( strcmp(chTime,"LAST") != 0 ) FeIIPunPop(save.ipPnunit[ipPun], true, (long int)save.punarg[ipPun][0], (long int)save.punarg[ipPun][1], lgFlag); } /* save spectra in fits format */ else if( strcmp(save.chSave[ipPun],"FITS") == 0 ) { if( strcmp(chTime,"LAST") == 0 ) { saveFITSfile( save.ipPnunit[ipPun], NUM_OUTPUT_TYPES ); } } /* save gammas (but without element) */ else if( strcmp(save.chSave[ipPun],"GAMt") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) { long ns; /* save photoionization rates, with the PUNCH GAMMAS command */ for( nelem=0; nelem < LIMELM; nelem++ ) { for( ion=0; ion <= nelem; ion++ ) { for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ ) { fprintf( save.ipPnunit[ipPun], "%3ld%3ld%3ld%10.2e%10.2e%10.2e", nelem+1, ion+1, ns+1, ionbal.PhotoRate_Shell[nelem][ion][ns][0], ionbal.PhotoRate_Shell[nelem][ion][ns][1] , ionbal.PhotoRate_Shell[nelem][ion][ns][2] ); for( j=0; j < t_yield::Inst().nelec_eject(nelem,ion,ns); j++ ) { fprintf( save.ipPnunit[ipPun], "%5.2f", t_yield::Inst().elec_eject_frac(nelem,ion,ns,j) ); } fprintf( save.ipPnunit[ipPun], "\n" ); } } } } } /* save gammas element, ion */ else if( strcmp(save.chSave[ipPun],"GAMe") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) { int ns; nelem = (long)save.punarg[ipPun][0]; ion = (long)save.punarg[ipPun][1]; /* valence shell */ ns = Heavy.nsShells[nelem][ion]-1; /* show what some of the ionization sources are */ GammaPrt( opac.ipElement[nelem][ion][ns][0] , opac.ipElement[nelem][ion][ns][1] , opac.ipElement[nelem][ion][ns][2] , save.ipPnunit[ipPun], ionbal.PhotoRate_Shell[nelem][ion][ns][0] , ionbal.PhotoRate_Shell[nelem][ion][ns][0]*0.1 ); } } else if( strcmp(save.chSave[ipPun],"GAUN") == 0 ) { /* save gaunt factors */ if( strcmp(chTime,"LAST") != 0 ) SaveGaunts(save.ipPnunit[ipPun]); } // generating the SAVE GRID output has been moved to cdPrepareExit() // to make sure that the output always records any type of failure } else { //no hit this branch, key should be in next lgNoHitFirstBranch = true; } // hack needed for code to compile with Visual Studio // keep this identical to the if-statement further up!! if( iterations.lgLastIt || !save.lgPunLstIter[ipPun] || ( lgAbort && strcmp(chTime,"LAST") == 0 ) || ( dynamics.lgTimeDependentStatic && dynamics.lgStatic_completed ) ) { if( strcmp(save.chSave[ipPun],"HISp") == 0 ) { /* save pressure history of current zone */ if( strcmp(chTime,"LAST") != 0 ) { /* note if pressure convergence failure occurred in history that follows */ if( !conv.lgConvPres ) { fprintf( save.ipPnunit[ipPun], "#PROBLEM Pressure not converged iter %li zone %li density-pressure follows:\n", iteration , nzone ); } /* note if temperature convergence failure occurred in history that follows */ if( !conv.lgConvTemp ) { fprintf( save.ipPnunit[ipPun], "#PROBLEM Temperature not converged iter %li zone %li density-pressure follows:\n", iteration , nzone ); } for( unsigned long k=0; k < conv.hist_pres_density.size(); ++k ) { /* save history of density - pressure, with correct pressure */ fprintf( save.ipPnunit[ipPun] , "%2li %4li\t%.5e\t%.5e\t%.5e\n", iteration, nzone, conv.hist_pres_density[k], conv.hist_pres_current[k], conv.hist_pres_error[k]); } } } else if( strcmp(save.chSave[ipPun],"HISt") == 0 ) { /* save temperature history of current zone */ if( strcmp(chTime,"LAST") != 0 ) { /* note if pressure convergence failure occurred in history that follows */ if( !conv.lgConvPres ) { fprintf( save.ipPnunit[ipPun], "#PROBLEM Pressure not converged iter %li zone %li temp heat cool follows:\n", iteration , nzone ); } /* note if temperature convergence failure occurred in history that follows */ if( !conv.lgConvTemp ) { fprintf( save.ipPnunit[ipPun], "#PROBLEM Temperature not converged iter %li zone %li temp heat cool follows:\n", iteration , nzone ); } for( unsigned long k=0; k < conv.hist_temp_temp.size(); ++k ) { /* save history of density - pressure, with correct pressure */ fprintf( save.ipPnunit[ipPun] , "%2li %4li\t%.5e\t%.5e\t%.5e\n", iteration, nzone, conv.hist_temp_temp[k], conv.hist_temp_heat[k], conv.hist_temp_cool[k]); } } } else if( strncmp(save.chSave[ipPun],"H2",2) == 0 ) { /* all save info on large H2 molecule include H2 PDR pdr */ save.whichDiatomToPrint[ipPun]->H2_PunchDo( save.ipPnunit[ipPun] , save.chSave[ipPun] , chTime, ipPun ); } else if( strcmp(save.chSave[ipPun],"HEAT") == 0 ) { /* save heating, routine in file of same name */ if( strcmp(chTime,"LAST") != 0 ) SaveHeat(save.ipPnunit[ipPun]); } else if( strncmp(save.chSave[ipPun],"HE",2) == 0 ) { /* various save helium commands */ /* save helium line wavelengths */ if( strcmp(save.chSave[ipPun] , "HELW") == 0 ) { if( strcmp(chTime,"LAST") == 0 ) { /* save helium & he-like wavelengths, first header */ fprintf( save.ipPnunit[ipPun], "Z\tElem\t2 1P->1 1S\t2 3P1->1 1S\t2 3P2->1 1S" "\t2 3S->1 1S\t2 3P2->2 3S\t2 3P1->2 3S\t2 3P0->2 3S" ); fprintf( save.ipPnunit[ipPun], "\n" ); for( nelem=ipHELIUM; nelemden/dense.gas_phase[ipHYDROGEN], findspecieslocal("H3+")->den/dense.gas_phase[ipHYDROGEN], findspecieslocal("H-")->den/dense.gas_phase[ipHYDROGEN] ); } } else if( strcmp(save.chSave[ipPun],"HYDi") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) { /* save hydrogen ionization * this will be total decays to ground */ RateInter = 0.; stage = iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ColIoniz*dense.eden*iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop(); fref = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc*iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop(); fout = iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop(); /* 06 aug 28, from numLevels_max to _local. */ for( ion=ipH2s; ion < iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_local; ion++ ) { /* this is total decays to ground */ RateInter += iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Aul()* (iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Pesc() + iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Pelec_esc() + iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Pdest()); /* total photo from all levels */ fref += iso_sp[ipH_LIKE][ipHYDROGEN].fb[ion].gamnc*iso_sp[ipH_LIKE][ipHYDROGEN].st[ion].Pop(); /* total col ion from all levels */ stage += iso_sp[ipH_LIKE][ipHYDROGEN].fb[ion].ColIoniz*dense.eden* iso_sp[ipH_LIKE][ipHYDROGEN].st[ion].Pop(); fout += iso_sp[ipH_LIKE][ipHYDROGEN].st[ion].Pop(); } /* make these relative to parent ion */ stage /= dense.xIonDense[ipHYDROGEN][1]; fref /= dense.xIonDense[ipHYDROGEN][1]; fout /= dense.xIonDense[ipHYDROGEN][1]; fprintf( save.ipPnunit[ipPun], "hion\t%4ld\t%.2e\t%.2e\t%.2e", nzone, /* photo and collision ion rates have units s-1 */ iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc, iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ColIoniz* dense.EdenHCorr, ionbal.RateRecomTot[ipHYDROGEN][0] ); fprintf( save.ipPnunit[ipPun], "\t%.2e", iso_sp[ipH_LIKE][ipHYDROGEN].RadRec_caseB ); fprintf( save.ipPnunit[ipPun], "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0], iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc/(ionbal.RateRecomTot[ipHYDROGEN][0]), iso_sp[ipH_LIKE][ipHYDROGEN].fb[1].RadRecomb[ipRecEsc], RateInter, fref/MAX2(1e-37,fout), stage/MAX2(1e-37,fout), /* simple H+ */ safe_div( iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc*dense.xIonDense[ipHYDROGEN][0], dense.eden*dense.xIonDense[ipHYDROGEN][1] ), secondaries.csupra[ipHYDROGEN][0]); GammaPrt(iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon,rfield.nflux,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipOpac, save.ipPnunit[ipPun],iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc* 0.05); } } else if( strcmp(save.chSave[ipPun],"HYDp") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) { /* save hydrogen populations * first give total atom and ion density [cm-3]*/ fprintf( save.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e", radius.depth_mid_zone, dense.xIonDense[ipHYDROGEN][0], dense.xIonDense[ipHYDROGEN][1] ); /* next give state-specific densities [cm-3] */ for( j=ipH1s; j < iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_local-1; j++ ) { fprintf( save.ipPnunit[ipPun], "\t%.2e", iso_sp[ipH_LIKE][ipHYDROGEN].st[j].Pop() ); } fprintf( save.ipPnunit[ipPun], "\n" ); } } else if( strcmp(save.chSave[ipPun],"HYDl") == 0 ) { if( strcmp(chTime,"LAST") == 0 ) { /* save hydrogen line * gives intensities and optical depths */ for( ipHi=1; ipHi>chng 04 oct 15, from nelem+2 to nelem+1 - array over read - * caught by PnH */ for( ion=0; ion= Heavy.nsShells[nelem][ion] ) break; /* array elements are shell, numb of electrons, element, 0 */ energy = t_ADfA::Inst().ph1(ips,nelem-ion,nelem,0); /* now print threshold with correct format */ if( energy < 10. ) { fprintf( save.ipPnunit[ipPun], "\t%6.3f", energy ); } else if( energy < 100. ) { fprintf( save.ipPnunit[ipPun], "\t%6.2f", energy ); } else if( energy < 1000. ) { fprintf( save.ipPnunit[ipPun], "\t%6.1f", energy ); } else { fprintf( save.ipPnunit[ipPun], "\t%6ld", (long)(energy) ); } } /* put cs at end of long line */ fprintf( save.ipPnunit[ipPun], "\n" ); } } } } } else if( strcmp(save.chSave[ipPun],"LINC") == 0 ) { /* save line cumulative */ if( strcmp(chTime,"LAST") != 0 ) { save_line(save.ipPnunit[ipPun],"PUNC", save.lgEmergent[ipPun]); } } else if( strcmp(save.chSave[ipPun],"LIND") == 0 ) { /* save line data, then stop */ SaveLineData(save.ipPnunit[ipPun]); } else if( strcmp(save.chSave[ipPun],"LINL") == 0 ) { /* save line labels */ bool lgPrintAll=false; /* LONG keyword on save line labels command sets this to 1 */ if( save.punarg[ipPun][0]>0. ) lgPrintAll = true; prt_LineLabels(save.ipPnunit[ipPun] , lgPrintAll ); } else if( strcmp(save.chSave[ipPun],"LINO") == 0 ) { if( strcmp(chTime,"LAST") == 0 ) { /* save line optical depths, routine is below, file static */ SaveLineStuff(save.ipPnunit[ipPun],"optical" , save.punarg[ipPun][0]); } } else if( strcmp(save.chSave[ipPun],"LINP") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) { static bool lgFirst=true; /* save line populations, need to do this twice if very first * call since first call to SaveLineStuff generates atomic parameters * rather than level pops, routine is below, file static */ SaveLineStuff(save.ipPnunit[ipPun],"populat" , save.punarg[ipPun][0]); if( lgFirst ) { lgFirst = false; SaveLineStuff(save.ipPnunit[ipPun],"populat" , save.punarg[ipPun][0]); } } } else if( strcmp(save.chSave[ipPun],"LINS") == 0 ) { /* save line emissivity */ if( strcmp(chTime,"LAST") != 0 ) { save_line(save.ipPnunit[ipPun],"PUNS", save.lgEmergent[ipPun]); } } else if( strcmp(save.chSave[ipPun],"LINR") == 0 ) { /* save line RT */ if( strcmp(chTime,"LAST") != 0 ) Save_Line_RT( save.ipPnunit[ipPun]); } else if( strcmp(save.chSave[ipPun],"LINA") == 0 ) { /* save line array */ if( strcmp(chTime,"LAST") == 0 ) { /* save out all lines with energies */ for( j=0; j < LineSave.nsum; j++ ) { if( LineSv[j].wavelength > 0. && LineSv[j].SumLine[0] > 0. ) { /* line energy, in units set with units option */ fprintf( save.ipPnunit[ipPun], "%12.5e", AnuUnit((realnum)RYDLAM/LineSv[j].wavelength) ); /* line label */ fprintf( save.ipPnunit[ipPun], "\t%4.4s ", LineSv[j].chALab ); /* wavelength */ prt_wl( save.ipPnunit[ipPun], LineSv[j].wavelength ); /* intrinsic intensity */ fprintf( save.ipPnunit[ipPun], "\t%8.3f", log10(SDIV(LineSv[j].SumLine[0]) ) + radius.Conv2PrtInten ); /* emergent line intensity, r recombination */ fprintf( save.ipPnunit[ipPun], "\t%8.3f", log10(SDIV(LineSv[j].SumLine[1]) ) + radius.Conv2PrtInten ); /* type of line, i for info, etc */ fprintf( save.ipPnunit[ipPun], " \t%c\n", LineSv[j].chSumTyp); } } } } else if( strcmp(save.chSave[ipPun],"LINI") == 0 ) { if( strcmp(chTime,"LAST") == 0 && (nzone/save.LinEvery)*save.LinEvery != nzone ) { /* this is the last zone * save line intensities - but do not do last zone twice */ SaveLineIntensity(save.ipPnunit[ipPun] , ipPun , save.punarg[ipPun][0] ); } else if( strcmp(chTime,"LAST") != 0 ) { /* following so we only save first zone if LinEvery reset */ if( (save.lgLinEvery && nzone == 1) || (nzone/save.LinEvery)*save.LinEvery == nzone ) { /* this is middle of calculation * save line intensities */ SaveLineIntensity(save.ipPnunit[ipPun] , ipPun , save.punarg[ipPun][0]); } } } else if( strcmp( save.chSave[ipPun],"LEIL") == 0) { /* some line intensities for the Leiden PDR, * but only do this when calculation is complete */ if( strcmp(chTime,"LAST") == 0 ) { double absval , rel; long int n; /* the lines we will find, * for a sample list of PDR lines look at LineList_PDR_H2.dat * in the cloudy data dir */ /* the number of H2 lines */ const int NLINE_H2 = 31; /* the number of lines which are not H2 */ const int NLINE_NOTH_H2 = 5; /* the labels and wavelengths for the lines that are not H2 */ char chLabel[NLINE_NOTH_H2][5]= {"C 2", "O 1", "O 1","C 1", "C 1" }; double Wl[NLINE_NOTH_H2]= {157.6 , 63.17 , 145.5 ,609.2 , 369.7 , }; /* these are wavelengths in microns, conv to Angstroms before call */ /* >>chng 05 sep 06, many of following wavelengths updated to agree * with output - apparently not updated when energies changed */ double Wl_H2[NLINE_H2]= {2.121 , 28.213, 17.03 , 12.28 , 9.662 , 8.024 , 6.907 , 6.107 , 5.510 , 5.051 , 4.693 , 4.408 , 4.180 , 3.996 , 3.845 , 3.723 , 3.625 , 3.546 , 3.485 , 3.437 , 3.403 , 3.380 , 3.368 , 3.365 , 3.371 , 3.387 , 3.410 , 3.441 , 3.485 , 3.542 , 3.603}; /* print a header for the lines */ for( n=0; nden, findspecieslocal("O2")->den, findspecieslocal("CH")->den, findspecieslocal("OH")->den, dense.eden, dense.xIonDense[ipHELIUM][1], dense.xIonDense[ipHYDROGEN][1], findspecieslocal("H3+")->den, colden.colden[ipCOL_H0], colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s], col_ci, col_cii, col_oi, findspecieslocal("CO")->column, findspecieslocal("O2")->column, findspecieslocal("CH")->column, findspecieslocal("OH")->column, colden.colden[ipCOL_elec], col_heii, colden.colden[ipCOL_Hp], colden.colden[ipCOL_H3p], hmi.H2_Solomon_dissoc_rate_used_H2g , gv.rate_h2_form_grains_used_total, hmi.H2_photodissoc_used_H2g, hmi.UV_Cont_rel2_Draine_DB96_depth, /* CO and C dissociation rate */ mole.findrk("PHOTON,CO=>C,O"), /* total CI ionization rate */ ionbal.PhotoRate_Shell[ipCARBON][0][2][0], /* total heating, erg cm-3 s-1 */ thermal.htot, /* total cooling, erg cm-3 s-1 */ thermal.ctot, /* GrnP grain photo heating */ thermal.heating[0][13], /* grain collisional cooling */ MAX2(0.,gv.GasCoolColl), /* grain collisional heating */ -1.*MIN2(0.,gv.GasCoolColl), /* COds - CO dissociation heating */ thermal.heating[0][9], /* H2dH-Heating due to H2 dissociation */ hmi.HeatH2Dish_used, /* H2vH-Heating due to collisions with H2 */ hmi.HeatH2Dexc_used , /* ChaT - charge transfer heating */ thermal.heating[0][24] , /* cosmic ray heating */ thermal.heating[1][6] , /* heating due to atoms of various heavy elements */ thermal.heating[ipMAGNESIUM][0], thermal.heating[ipSULPHUR][0], thermal.heating[ipSILICON][0], thermal.heating[ipIRON][0], thermal.heating[ipSODIUM][0], thermal.heating[ipALUMINIUM][0], thermal.heating[ipCARBON][0], TauLines[ipT610].Coll().cool(), TauLines[ipT370].Coll().cool(), TauLines[ipT157].Coll().cool(), TauLines[ipT63].Coll().cool(), TauLines[ipT146].Coll().cool() ); } } # ifdef USE_NLTE7 else if( strcmp( save.chSave[ipPun],"NLTE") == 0) { if( strcmp(chTime,"LAST") == 0 ) { //Generate output for NLTE7 conference runNLTE(ipPun); } } # endif /* Print out the FeII energy level file into Stout format*/ else if( strcmp( save.chSave[ipPun],"LY1") == 0) { static bool runonce = true; if (runonce ) { for( long ipHi=0; ipHi < FeII.nFeIILevel_malloc; ipHi++ ) { ASSERT(FeII.FeIINRGs[ipHi] >= 0.); ASSERT(FeII.FeIISTWT[ipHi] >= 0.); fprintf(save.ipPnunit[ipPun], "%li\t%10.3f\t%3.1f\n",ipHi+1,FeII.FeIINRGs[ipHi],FeII.FeIISTWT[ipHi]); } runonce = false; fprintf(save.ipPnunit[ipPun], "-1\n"); } } /* Print out the FeII transition probablility file into Stout format*/ else if( strcmp( save.chSave[ipPun],"LY2") == 0) { static bool runonce = true; if (runonce ) { for( long ipLo = 0; ipLo < FeII.nFeIILevel_malloc; ipLo++) { for( long ipHi=0; ipHi < FeII.nFeIILevel_malloc; ipHi++ ) { if( FeII.FeIIAul[ipHi][ipLo] >= 0. && FeII.FeIIAul[ipHi][ipLo] != 1e-5f && FeII.FeIIAul[ipHi][ipLo] != 1e-6f ) { fprintf(save.ipPnunit[ipPun], "A\t%li\t%li\t%8.2e\n", ipLo+1,ipHi+1,FeII.FeIIAul[ipHi][ipLo]); } } } runonce = false; fprintf(save.ipPnunit[ipPun], "-1\n"); } } /* Print out the FeII Collision data file into Stout format*/ else if( strcmp( save.chSave[ipPun],"LY3") == 0) { static realnum tt[8]={1e3f,3e3f,5e3f,7e3f,1e4f,12e3f,15e3f,2e4f}; static bool runonce = true; if (runonce ) { fprintf(save.ipPnunit[ipPun], "TEMP\t8"); for( int k = 0; k < 8; k++) { fprintf(save.ipPnunit[ipPun], "\t%8.2e",tt[k]); } fprintf(save.ipPnunit[ipPun], "\n"); for( long ipHi = 1; ipHi < FeII.nFeIILevel_malloc; ipHi++) { for( long ipLo=0; ipLo < ipHi; ipLo++ ) { bool skipLine = false; for( int k = 0; k < 8; k++) { if( FeII.FeIIColl[ipHi][ipLo][k] == -2. || FeII.FeIIColl[ipHi][ipLo][k] == 0.) { skipLine = true; break; } } if( !skipLine ) { fprintf(save.ipPnunit[ipPun], "CSELECTRON\t%li\t%li",ipLo+1,ipHi+1); for( int k = 0; k < 8; k++) { fprintf(save.ipPnunit[ipPun], "\t%8.2e",FeII.FeIIColl[ipHi][ipLo][k]); } fprintf(save.ipPnunit[ipPun], "\n"); } } } runonce = false; fprintf(save.ipPnunit[ipPun], "-1\n"); } } else if( strcmp( save.chSave[ipPun],"LLST") == 0) { /* save linelist command - do on last iteration */ if( strcmp(chTime,"LAST") == 0 ) { fprintf( save.ipPnunit[ipPun], "iteration %li" , iteration ); if( save.punarg[ipPun][1] )// column print fprintf( save.ipPnunit[ipPun], "\n" ); /* -1 is flag saying that this save command was not set */ if( save.nLineList[ipPun] < 0 ) TotalInsanity(); int LineType = 0; if( save.lgEmergent[ipPun] ) LineType = 1; if( save.lgCumulative[ipPun] ) LineType += 2; /* loop over all lines in the file we read */ for( j=0; j 0 ) PrtQuantity = pow(10. , absolute); else PrtQuantity = relative; // column mode, print label if( save.punarg[ipPun][1] ) { /* if taking ratio then put div sign between pairs */ if( save.lgLineListRatio[ipPun] && is_odd(j) ) fprintf( save.ipPnunit[ipPun] , "/" ); fprintf( save.ipPnunit[ipPun], "%s ", save.chLineListLabel[ipPun][j] ); char chTemp[MAX_HEADER_SIZE]; sprt_wl( chTemp, save.wlLineList[ipPun][j] ); fprintf( save.ipPnunit[ipPun], "%s ", chTemp ); } /* if taking ratio print every other line as ratio * with previous line */ if( save.lgLineListRatio[ipPun] ) { /* do line pair ratios */ static double SaveQuantity = 0; if( is_odd(j) ) fprintf( save.ipPnunit[ipPun], "\t%.4e" , SaveQuantity / SDIV( PrtQuantity ) ); else SaveQuantity = PrtQuantity; } else { fprintf( save.ipPnunit[ipPun], "\t%.4e" , PrtQuantity ); } // column printout, but check if first of pair if( save.punarg[ipPun][1] ) { if( !save.lgLineListRatio[ipPun] || is_odd(j) ) fprintf( save.ipPnunit[ipPun], "\n" ); } } fprintf( save.ipPnunit[ipPun], "\n" ); } } else if( strcmp( save.chSave[ipPun],"CHRT") == 0) { /* save chemistry rates command */ if( strcmp(chTime,"LAST") != 0 ) { bool lgHeader, lgData; if( save.lgPunHeader[ipPun] ) { lgHeader = true; lgData = false; mole_punch(save.ipPnunit[ipPun],save.optname[ipPun].c_str(),save.chSaveArgs[ipPun],lgHeader,lgData,radius.depth_mid_zone); } save.lgPunHeader[ipPun] = false; lgHeader = false; lgData = true; mole_punch(save.ipPnunit[ipPun],save.optname[ipPun].c_str(),save.chSaveArgs[ipPun],lgHeader,lgData,radius.depth_mid_zone); } } else if( strcmp(save.chSave[ipPun],"MAP ") == 0 ) { /* do the map now if we are at the zone, or if this * is the LAST call to this routine and map not done yet */ if( !hcmap.lgMapDone && (nzone == hcmap.MapZone || strcmp(chTime,"LAST") == 0 ) ) { lgTlkSav = called.lgTalk; called.lgTalk = cpu.i().lgMPI_talk(); hcmap.lgMapBeingDone = true; map_do(save.ipPnunit[ipPun] , " map"); called.lgTalk = lgTlkSav; } } else if( strcmp(save.chSave[ipPun],"MOLE") == 0 ) { if( save.lgPunHeader[ipPun] ) { fprintf( save.ipPnunit[ipPun], "#molecular species will follow:\n"); fprintf( save.ipPnunit[ipPun], "#depth\tAV(point)\tAV(extend)\tCO diss rate\tC recom rate"); for(i=0; ilabel.c_str() ); } fprintf ( save.ipPnunit[ipPun], "\n"); save.lgPunHeader[ipPun] = false; } if( strcmp(chTime,"LAST") != 0 ) { /* molecules, especially for PDR, first give radius */ fprintf( save.ipPnunit[ipPun], "%.5e\t" , radius.depth_mid_zone ); /* visual extinction of point source (star)*/ fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_point); /* visual extinction of an extended source (like a PDR)*/ fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended); /* carbon monoxide photodissociation rate */ fprintf( save.ipPnunit[ipPun], "%.5e\t" , mole.findrk("PHOTON,CO=>C,O") ); /* carbon recombination rate */ fprintf( save.ipPnunit[ipPun], "%.5e" , ionbal.RateRecomTot[ipCARBON][0] ); /* now do all the molecules */ for(j=0; j 0. ) /* possible lower bounds to energy range - will be zero if not set */ j = ipFineCont( save.punarg[ipPun][0] ); else j = 0; /* upper limit */ if( save.punarg[ipPun][1]> 0. ) nu_hi = ipFineCont( save.punarg[ipPun][1]); else nu_hi = rfield.nfine; /* we will bring nskip cells together into one printed * number to make output smaller - default is 10 */ nskip = (long)save.punarg[ipPun][2]; do { realnum sum1 = rfield.fine_opt_depth[j]; realnum sum2 = rfield.fine_opac_zone[j]; /* want to report the central wavelength of the cell */ realnum xnu = rfield.fine_anu[j]; for( jj=1; jj 0. ) fprintf( save.ipPnunit[ipPun], "%12.6e\t%.3e\t%.3e\n", AnuUnit(xnu/nskip), sum1/nskip , sum2/nskip); j += nskip; }while( j < nu_hi ); } } else if( strcmp(save.chSave[ipPun]," OTS") == 0 ) { ConMax = 0.; xLinMax = 0.; opConSum = 0.; opLinSum = 0.; ipLinMax = 1; ipConMax = 1; for( j=0; j < rfield.nflux; j++ ) { opConSum += rfield.otscon[j]*opac.opacity_abs[j]; opLinSum += rfield.otslin[j]*opac.opacity_abs[j]; if( rfield.otslin[j]*opac.opacity_abs[j] > xLinMax ) { xLinMax = rfield.otslin[j]*opac.opacity_abs[j]; ipLinMax = j+1; } if( rfield.otscon[j]*opac.opacity_abs[j] > ConMax ) { ConMax = rfield.otscon[j]*opac.opacity_abs[j]; ipConMax = j+1; } } fprintf( save.ipPnunit[ipPun], "tot con lin=%.2e%.2e lin=%.4s%.4e%.2e con=%.4s%.4e%.2e\n", opConSum, opLinSum, rfield.chLineLabel[ipLinMax-1] , rfield.anu[ipLinMax-1], xLinMax, rfield.chContLabel[ipConMax-1] , rfield.anu[ipConMax-1], ConMax ); } else if( strcmp(save.chSave[ipPun],"OVER") == 0 ) { /* save overview * this is the floor for the smallest ionization fractions printed */ double toosmall = SMALLFLOAT , hold; /* overview of model results, * depth, te, hden, eden, ion fractions H, He, c, O */ if( strcmp(chTime,"LAST") != 0 ) { /* print the depth */ fprintf( save.ipPnunit[ipPun], "%.5e\t", radius.depth_mid_zone ); /* temperature, heating */ if(dynamics.Cool() > dynamics.Heat()) { fprintf( save.ipPnunit[ipPun], "%.4f\t%.3f", log10(phycon.te), log10(SDIV(thermal.htot-dynamics.Heat())) ); } else { double diff = fabs(thermal.htot-dynamics.Cool()); fprintf( save.ipPnunit[ipPun], "%.4f\t%.3f", log10(phycon.te), log10(SDIV(diff)) ); } /* hydrogen and electron densities */ fprintf( save.ipPnunit[ipPun], "\t%.4f\t%.4f", log10(dense.gas_phase[ipHYDROGEN]), log10(dense.eden) ); /* molecular fraction of hydrogen */ fprintf( save.ipPnunit[ipPun], "\t%.4f", /*log10(MAX2(toosmall,2.*findspecies("H2")->den/dense.gas_phase[ipHYDROGEN])) );*/ log10(MAX2(toosmall,2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN])) ); /* ionization fractions of hydrogen */ fprintf( save.ipPnunit[ipPun], "\t%.4f\t%.4f", log10(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN])), log10(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN])) ); /* ionization fractions of helium */ for( j=1; j <= 3; j++ ) { double arg1 = SDIV(dense.gas_phase[ipHELIUM]); arg1 = MAX2(toosmall,dense.xIonDense[ipHELIUM][j-1]/arg1 ); fprintf( save.ipPnunit[ipPun], "\t%.4f", log10(arg1) ); } /* carbon monoxide molecular fraction of CO */ hold = SDIV(dense.gas_phase[ipCARBON]); hold = findspecieslocal("CO")->den/hold; hold = MAX2(toosmall, hold ); fprintf( save.ipPnunit[ipPun], "\t%.4f", log10(hold) ); /* ionization fractions of carbon */ for( j=1; j <= 4; j++ ) { hold = SDIV(dense.gas_phase[ipCARBON]); hold = MAX2(toosmall,dense.xIonDense[ipCARBON][j-1]/hold); fprintf( save.ipPnunit[ipPun], "\t%.4f", log10(hold) ); } /* ionization fractions of oxygen */ for( j=1; j <= 6; j++ ) { hold = SDIV(dense.gas_phase[ipOXYGEN]); hold = MAX2(toosmall,dense.xIonDense[ipOXYGEN][j-1]/hold); fprintf( save.ipPnunit[ipPun], "\t%.4f", log10(hold) ); } // molecular fraction of H2O hold = SDIV(dense.gas_phase[ipOXYGEN]); hold = findspecieslocal("H2O")->den/hold; hold = MAX2(toosmall, hold ); fprintf( save.ipPnunit[ipPun], "\t%.4f", log10(hold) ); /* visual extinction of point source (star)*/ fprintf( save.ipPnunit[ipPun], "\t%.2e" , rfield.extin_mag_V_point); /* visual extinction of an extended source (like a PDR)*/ fprintf( save.ipPnunit[ipPun], "\t%.2e\n" , rfield.extin_mag_V_extended); } } else if( strcmp(save.chSave[ipPun]," PDR") == 0 ) { /* this is the save PDR command */ if( strcmp(chTime,"LAST") != 0 ) { /* convert optical depth at wavelength of V filter * into magnitudes of extinction */ /* >>chyng 03 feb 25, report extinction to illuminated face, * rather than total extinction which included far side when * sphere was set */ /*av = opac.TauTotalGeo[0][rfield.ipV_filter-1]*1.08574;*/ fprintf( save.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t", radius.depth_mid_zone, /* total hydrogen column density, all forms */ colden.colden[ipCOL_HTOT], phycon.te, /* fraction of H that is atomic */ dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN], /* ratio of n(H2) to total H, == 0.5 when fully molecular */ 2.*findspecieslocal("H2")->den/dense.gas_phase[ipHYDROGEN], 2.*findspecieslocal("H2*")->den/dense.gas_phase[ipHYDROGEN], /* atomic to total carbon */ dense.xIonDense[ipCARBON][0]/SDIV(dense.gas_phase[ipCARBON]), findspecieslocal("CO")->den/SDIV(dense.gas_phase[ipCARBON]), findspecieslocal("H2O")->den/SDIV(dense.gas_phase[ipOXYGEN]), /* hmi.UV_Cont_rel2_Habing_TH85 is field relative to Habing background, dimensionless */ hmi.UV_Cont_rel2_Habing_TH85_depth); /* visual extinction due to dust alone, of point source (star)*/ fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_point); /* visual extinction due to dust alone, of an extended source (like a PDR)*/ fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended); /* visual extinction (all sources) of a point source (like a PDR)*/ fprintf( save.ipPnunit[ipPun], "%.2e\n" , opac.TauAbsGeo[0][rfield.ipV_filter] ); } } /* performance characteristics per zone */ else if( strcmp(save.chSave[ipPun],"PERF") == 0 ) { if( strcmp(chTime,"LAST") != 0 ) { static double ElapsedTime , ZoneTime; if( nzone<=1 ) { ElapsedTime = cdExecTime(); ZoneTime = 0.; } else { double t = cdExecTime(); ZoneTime = t - ElapsedTime; ElapsedTime = t; } /* zone, time for this zone, elapsed time */ fprintf( save.ipPnunit[ipPun], " %ld\t%.3f\t%.2f\t%li", nzone, ZoneTime , ElapsedTime, conv.nPres2Ioniz ); // print various loop counters for( long i=0; i>chng 06 apr 19, subtract pinzon the acceleration added in this zone * since is not total at outer edge of zone, above is at inner edge */ pressure.PresTotlInit + pressure.PresInteg - pressure.pinzon, /*E 5 pgas (0) */ pressure.PresTotlInit, /*F 6 Pgas */ pressure.PresGasCurr, /*G 7 Pram */ pressure.PresRamCurr, /*H 8 P rad in lines */ pressure.pres_radiation_lines_curr, /*I 9 Pinteg subtract continuum rad pres which has already been added on */ pressure.PresInteg - pressure.pinzon, /*J 10 V(wind km/s) wind speed in km/s */ wind.windv/1e5, /*K cad(km/s) sound speed in km/s */ timesc.sound_speed_adiabatic/1e5, /* the magnetic pressure */ magnetic.pressure , /* the local turbulent velocity in km/s */ DoppVel.TurbVel/1e5 , /* turbulent pressure */ pressure.PresTurbCurr*DoppVel.lgTurb_pressure , /* gravitational pressure */ pressure.IntegRhoGravity, // the integral of electron scattering acceleration in // the absence of any absorptio, minus acceleration in current // zone, which has been added in - done this way, result is // zero in first zonen pressure.PresIntegElecThin-pressure.pinzon_PresIntegElecThin, // is this converged? TorF(conv.lgConvPres) ); } } else if( strcmp(save.chSave[ipPun],"PREL") == 0 ) { /* line pressure contributors */ fprintf( save.ipPnunit[ipPun], "%.5e\t%.3e\t%.3e\t", /*A 1 #P depth */ radius.depth_mid_zone , pressure.PresTotlCurr, pressure.pres_radiation_lines_curr/SDIV(pressure.PresTotlCurr) ); PrtLinePres(save.ipPnunit[ipPun]); } else if( save.chSave[ipPun][0]=='R' ) { /* work around internal limits to Microsoft vs compiler */ if( strcmp(save.chSave[ipPun],"RADI") == 0 ) { /* save radius information for all zones */ if( strcmp(chTime,"LAST") != 0 ) { fprintf( save.ipPnunit[ipPun], "%ld\t%.5e\t%.4e\t%.4e\n", nzone, radius.Radius_mid_zone, radius.depth_mid_zone, radius.drad ); } } else if( strcmp(save.chSave[ipPun],"RADO") == 0 ) { /* save radius information for only the last zone */ if( strcmp(chTime,"LAST") == 0 ) { fprintf( save.ipPnunit[ipPun], "%ld\t%.5e\t%.4e\t%.4e\n", nzone, radius.Radius_mid_zone, radius.depth_mid_zone, radius.drad ); } } else if( strcmp(save.chSave[ipPun],"RESU") == 0 ) { /* save results of the calculation */ if( strcmp(chTime,"LAST") == 0 ) SaveResults(save.ipPnunit[ipPun]); } else { /* this can't happen */ TotalInsanity(); } } else if( strcmp(save.chSave[ipPun],"SECO") == 0 ) { /* save secondary ionization */ if( strcmp(chTime,"LAST") != 0 ) fprintf(save.ipPnunit[ipPun], "%.5e\t%.3e\t%.3e\t%.3e\n", radius.depth , secondaries.csupra[ipHYDROGEN][0], secondaries.csupra[ipHYDROGEN][0]*2.02, secondaries.x12tot ); } else if( strcmp(save.chSave[ipPun],"SOUS") == 0 ) { /* full spectrum of continuum source function at 1 depth * command was "save source spectrum" */ if( strcmp(chTime,"LAST") != 0 ) { limit = MIN2(rfield.ipMaxBolt,rfield.nflux); for( j=0; j < limit; j++ ) { fprintf( save.ipPnunit[ipPun], "%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n", AnuUnit(rfield.AnuOrg[j]), rfield.ConEmitLocal[nzone][j]/rfield.widflx[j], opac.opacity_abs[j], rfield.ConSourceFcnLocal[nzone][j], rfield.ConSourceFcnLocal[nzone][j]/plankf(j) , safe_div(rfield.ConSourceFcnLocal[nzone][j],rfield.flux[0][j])); } } } else if( strcmp(save.chSave[ipPun],"SOUD") == 0 ) { /* parts of continuum source function vs depth * command was save source function depth */ j = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon + 2; fprintf( save.ipPnunit[ipPun], "%.4e\t%.4e\t%.4e\t%.4e\n", opac.TauAbsFace[j-1], rfield.ConEmitLocal[nzone][j-1]/rfield.widflx[j-1]/MAX2(1e-35,opac.opacity_abs[j-1]), rfield.otscon[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1], rfield.otscon[iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]/opac.opacity_abs[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1] ); } /* this is save special option */ else if( strcmp(save.chSave[ipPun],"SPEC") == 0 ) { SaveSpecial(save.ipPnunit[ipPun],chTime); } /* this is save species option */ else if( strcmp(save.chSave[ipPun],"SPCS") == 0 ) { if( ( strcmp(chTime,"LAST") != 0 && strcmp(save.chSaveArgs[ipPun],"COLU") != 0 ) || ( strcmp(chTime,"LAST") == 0 && strcmp(save.chSaveArgs[ipPun],"COLU") == 0 ) ) SaveSpecies(save.ipPnunit[ipPun] , ipPun); } else if( strcmp(save.chSave[ipPun],"TEMP") == 0 ) { static double deriv_old=-1; double deriv=-1. , deriv_sec; /* temperature and its derivatives */ fprintf( save.ipPnunit[ipPun], "%.5e\t%.4e\t%.2e", radius.depth_mid_zone, phycon.te, thermal.dCooldT ); /* if second zone then have one deriv */ if( nzone >1 ) { deriv = (phycon.te - struc.testr[nzone-2])/ radius.drad; fprintf( save.ipPnunit[ipPun], "\t%.2e", deriv ); /* if third zone then have second deriv */ if( nzone > 2 ) { deriv_sec = (deriv-deriv_old)/ radius.drad; fprintf( save.ipPnunit[ipPun], "\t%.2e", deriv_sec ); } deriv_old = deriv; } fprintf( save.ipPnunit[ipPun], "\n"); } /* time dependent model */ else if( strcmp(save.chSave[ipPun],"TIMD") == 0 ) { if( strcmp(chTime,"LAST") == 0 ) DynaPunchTimeDep( save.ipPnunit[ipPun] , "END" ); } /* execution time per zone */ else if( strcmp(save.chSave[ipPun],"XTIM") == 0 ) { static double ElapsedTime , ZoneTime; if( nzone<=1 ) { ElapsedTime = cdExecTime(); ZoneTime = 0.; } else { double t = cdExecTime(); ZoneTime = t - ElapsedTime; ElapsedTime = t; } /* zone, time for this zone, elapsed time */ fprintf( save.ipPnunit[ipPun], " %ld\t%.3f\t%.2f\n", nzone, ZoneTime , ElapsedTime ); } else if( strcmp(save.chSave[ipPun],"TPRE") == 0 ) { /* temperature and its predictors, turned on with save tprid */ fprintf( save.ipPnunit[ipPun], "%5ld %11.4e %11.4e %11.4e %g\n", nzone, phycon.TeInit, phycon.TeProp, phycon.te, (phycon.TeProp- phycon.te)/phycon.te ); } else if( strcmp(save.chSave[ipPun],"WIND") == 0 ) { /* wind velocity, radiative acceleration, and ratio total * to electron scattering acceleration */ /* first test only save last zone */ if( (save.punarg[ipPun][0] == 0 && strcmp(chTime,"LAST") == 0) || /* this test save all zones */ (save.punarg[ipPun][0] == 1 && strcmp(chTime,"LAST") != 0 ) ) { fprintf( save.ipPnunit[ipPun], "%.5e\t%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n", radius.Radius_mid_zone, radius.depth_mid_zone, wind.windv, wind.AccelTotalOutward, wind.AccelLine, wind.AccelCont , wind.fmul , wind.AccelGravity ); } } else if( strcmp(save.chSave[ipPun],"XATT") == 0 ) { /* attenuated incident continuum */ ASSERT( grid.lgOutputTypeOn[2] ); if( strcmp(chTime,"LAST") == 0 ) { if( grid.lgGrid ) saveFITSfile( save.ipPnunit[ipPun], 2 ); else { fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" ); cdEXIT(EXIT_FAILURE); } } } else if( strcmp(save.chSave[ipPun],"XRFI") == 0 ) { /* reflected incident continuum */ ASSERT( grid.lgOutputTypeOn[3] ); if( strcmp(chTime,"LAST") == 0 ) { if( grid.lgGrid ) saveFITSfile( save.ipPnunit[ipPun], 3 ); else { fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" ); cdEXIT(EXIT_FAILURE); } } } else if( strcmp(save.chSave[ipPun],"XINC") == 0 ) { /* incident continuum */ ASSERT( grid.lgOutputTypeOn[1] ); if( strcmp(chTime,"LAST") == 0 ) { if( grid.lgGrid ) saveFITSfile( save.ipPnunit[ipPun], 1 ); else { fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" ); cdEXIT(EXIT_FAILURE); } } } else if( strcmp(save.chSave[ipPun],"XDFR") == 0 ) { /* reflected diffuse continuous emission */ ASSERT( grid.lgOutputTypeOn[5] ); if( strcmp(chTime,"LAST") == 0 ) { if( grid.lgGrid ) saveFITSfile( save.ipPnunit[ipPun], 5 ); else { fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" ); cdEXIT(EXIT_FAILURE); } } } else if( strcmp(save.chSave[ipPun],"XDFO") == 0 ) { /* diffuse continuous emission outward */ ASSERT( grid.lgOutputTypeOn[4] ); if( strcmp(chTime,"LAST") == 0 ) { if( grid.lgGrid ) saveFITSfile( save.ipPnunit[ipPun], 4 ); else { fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" ); cdEXIT(EXIT_FAILURE); } } } else if( strcmp(save.chSave[ipPun],"XLNR") == 0 ) { /* reflected lines */ ASSERT( grid.lgOutputTypeOn[7] ); if( strcmp(chTime,"LAST") == 0 ) { if( grid.lgGrid ) saveFITSfile( save.ipPnunit[ipPun], 7 ); else { fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" ); cdEXIT(EXIT_FAILURE); } } } else if( strcmp(save.chSave[ipPun],"XLNO") == 0 ) { /* outward lines */ ASSERT( grid.lgOutputTypeOn[6] ); if( strcmp(chTime,"LAST") == 0 ) { if( grid.lgGrid ) saveFITSfile( save.ipPnunit[ipPun], 6 ); else { fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" ); cdEXIT(EXIT_FAILURE); } } } else if( strcmp(save.chSave[ipPun],"XREF") == 0 ) { /* total reflected, lines and continuum */ ASSERT( grid.lgOutputTypeOn[9] ); if( strcmp(chTime,"LAST") == 0 ) { if( grid.lgGrid ) saveFITSfile( save.ipPnunit[ipPun], 9 ); else { fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" ); cdEXIT(EXIT_FAILURE); } } } else if( strcmp(save.chSave[ipPun],"XTOT") == 0 ) { /* total spectrum, reflected plus transmitted */ ASSERT( grid.lgOutputTypeOn[0] ); if( strcmp(chTime,"LAST") == 0 ) { if( grid.lgGrid ) saveFITSfile( save.ipPnunit[ipPun], 0 ); else { fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" ); cdEXIT(EXIT_FAILURE); } } } else if( strcmp(save.chSave[ipPun],"XTRN") == 0 ) { /* total outward, lines and continuum */ ASSERT( grid.lgOutputTypeOn[8] ); if( strcmp(chTime,"LAST") == 0 ) { if( grid.lgGrid ) saveFITSfile( save.ipPnunit[ipPun], 8 ); else { fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" ); cdEXIT(EXIT_FAILURE); } } } else if( strcmp(save.chSave[ipPun],"XSPM") == 0 ) { /* exp(-tau) to the illuminated face */ ASSERT( grid.lgOutputTypeOn[10] ); if( strcmp(chTime,"LAST") == 0 ) { if( grid.lgGrid ) saveFITSfile( save.ipPnunit[ipPun], 10 ); else { fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" ); cdEXIT(EXIT_FAILURE); } } } // termination of second set of nested if's // error if we have not matched key /* there are a few "save" commands that are handled elsewhere * save dr is an example. These will have lgRealSave set false */ // lgNoHitFirstBranch says did not find in previous nest of if's else if( save.lgRealSave[ipPun] && lgNoHitFirstBranch ) { /* this is insanity, internal flag set in ParseSave not seen here */ fprintf( ioQQQ, " PROBLEM DISASTER SaveDo does not recognize flag %4.4s set by ParseSave. This is impossible.\n", save.chSave[ipPun] ); TotalInsanity(); } /* print special hash string to separate out various iterations * chTime is LAST on last iteration * save.lgHashEndIter flag is true by default, set false * with "no hash" keyword on save command * save.lg_separate_iterations is true by default, set false * when save time dependent calcs since do not want special * character between time steps * grid.lgGrid is only true when doing a grid of calculations */ if( strcmp(chTime,"LAST") == 0 && !(iterations.lgLastIt && !grid.lgGrid ) && save.lgHashEndIter[ipPun] && save.lg_separate_iterations[ipPun] && !save.lgFITS[ipPun] ) { if( dynamics.lgTimeDependentStatic && strcmp( save.chHashString , "TIME_DEP" )==0 ) { fprintf( save.ipPnunit[ipPun], "\"time=%f\n", dynamics.time_elapsed ); } else { fprintf( save.ipPnunit[ipPun], "%s", save.chHashString ); if( grid.lgGrid && ( iterations.lgLastIt || lgAbort ) ) fprintf( save.ipPnunit[ipPun], " GRID_DELIMIT -- grid%09ld", optimize.nOptimiz ); fprintf( save.ipPnunit[ipPun], "\n" ); } } if( save.lgFLUSH ) fflush( save.ipPnunit[ipPun] ); } } return; } /*SaveLineIntensity produce the 'save lines intensity' output */ STATIC void SaveLineIntensity(FILE * ioPUN, long int ipPun , realnum Threshold ) { long int i; DEBUG_ENTRY( "SaveLineIntensity()" ); /* used to save out all the emission line intensities * first initialize the line image reader */ fprintf( ioPUN, "**********************************************************************************************************************************\n" ); input.echo(ioPUN); /* now print any cautions or warnings */ cdWarnings( ioPUN); cdCautions( ioPUN); fprintf( ioPUN, "zone=%5ld\n", nzone ); fprintf( ioPUN, "**********************************************************************************************************************************\n" ); fprintf( ioPUN, "begin emission lines\n" ); /* only save non-zero intensities */ SaveResults1Line(ioPUN," ",0,0.,"Start"); // check whether intrinsic or emergent line emissivity bool lgEmergent = false; if( save.punarg[ipPun][0] > 0 ) lgEmergent = true; for( i=0; i < LineSave.nsum; i++ ) { // Threshold is zero by default on save line intensity, // all option sets to negative number so that we report all lines if( LineSv[i].SumLine[lgEmergent] > Threshold ) { SaveResults1Line(ioPUN,(char*)LineSv[i].chALab,LineSv[i].wavelength, LineSv[i].SumLine[save.lgEmergent[ipPun]], "Line "); } } SaveResults1Line(ioPUN," ",0,0.,"Flush"); fprintf( ioPUN, " \n" ); fprintf( ioPUN, "**********************************************************************************************************************************\n" ); return; } /* lgSaveOpticalDepths true says save optical depths */ static bool lgPopsFirstCall , lgSaveOpticalDepths; /*SaveLineStuff save optical depths or source functions for all transferred lines */ STATIC void SaveLineStuff( FILE * ioPUN, const char *chJob , realnum xLimit ) { long int nelem , ipLo , ipHi , i , ipISO; long index = 0; static bool lgFirst=true; DEBUG_ENTRY( "SaveLineStuff()" ); /*find out which job this is and set a flag to use later */ if( strcmp( &*chJob , "optical" ) == 0 ) { /* save line optical depths */ lgSaveOpticalDepths = true; lgPopsFirstCall = false; } else if( strcmp( &*chJob , "populat" ) == 0 ) { lgSaveOpticalDepths = false; /* level population information */ if( lgFirst ) { lgPopsFirstCall = true; fprintf(ioPUN,"index\tAn.ion\tgLo\tgUp\tE(wn)\tgf\n"); lgFirst = false; } else { lgPopsFirstCall = false; } } else { fprintf( ioQQQ, " insane job in SaveLineStuff =%s\n", &*chJob ); cdEXIT(EXIT_FAILURE); } /* loop over all lines, calling put1Line to create info (routine located below) */ /* hydrogen like lines */ /* >>chng 02 may 16, had been explicit H and He-like loops */ for( ipISO=ipH_LIKE; ipISO>chng 02 aug 23, for he-like, had starting on much too high a level since * index was number of levels - caught by Adrian Turner */ /* now output extra line lines, starting one above those already done above */ /*for( ipHi=iso_sp[ipISO][nelem].numLevels_max; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )*/ /* 06 aug 28, from numLevels_max to _local. */ for( ipHi=iso_sp[ipISO][nelem].st[iso_sp[ipISO][nelem].numLevels_local-1].n()+1; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ ) { ++index; Save1Line( ExtraLymanLines[ipISO][nelem][ipExtraLymanLines[ipISO][nelem][ipHi]], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[nelem]) ); } } } } } /* index from 1 to skip over dummy line */ for( i=1; i < nLevel1; i++ ) { ++index; Save1Line( TauLines[i], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[(*TauLines[i].Hi()).nelem()-1]) ); } for( i=0; i < nWindLine; i++ ) { if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO ) { ++index; Save1Line( TauLine2[i], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[(*TauLine2[i].Hi()).nelem()-1]) ); } } for( i=0; i < nUTA; i++ ) { ++index; Save1Line( UTALines[i], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[(*UTALines[i].Hi()).nelem()-1]) ); } /* do optical depths of FeII lines, if large atom is turned on */ FeIIPunchLineStuff( ioPUN , xLimit , index); /* save optical depths of H2 lines */ h2.H2_PunchLineStuff( ioPUN , xLimit , index); /*fprintf(ioPUN, "##################################\n"); */ fprintf( ioPUN , "%s\n",save.chHashString ); return; } /*Save1Line called by SaveLineStuff to produce output for one line */ void Save1Line( const TransitionProxy& t , FILE * ioPUN , realnum xLimit , long index, realnum DopplerWidth ) { if( lgSaveOpticalDepths ) { /* optical depths, no special first time, only print them */ if( t.Emis().TauIn() >= xLimit ) { /* label like "C 4" or "H 1" */ fprintf( ioPUN , "%2.2s%2.2s\t", elementnames.chElementSym[(*t.Hi()).nelem()-1] , elementnames.chIonStage[(*t.Hi()).IonStg()-1] ); /* print wavelengths, either line in main printout labels, * or in various units in exponential notation - prt_wl is in prt.c */ if( strcmp( save.chConPunEnr[save.ipConPun], "labl" )== 0 ) { prt_wl( ioPUN , t.WLAng() ); } else { /* this converts energy in Rydbergs into any of the other units */ fprintf( ioPUN , "%.7e", AnuUnit((realnum)(t.EnergyRyd())) ); } /* print the optical depth */ fprintf( ioPUN , "\t%.3f", t.Emis().TauIn() ); /* damping constant */ fprintf(ioPUN, "\t%.3e", t.Emis().dampXvel() / DopplerWidth ); fprintf(ioPUN, "\n"); } } else if( lgPopsFirstCall ) { char chLbl[11]; /* first call to line populations, print atomic parameters and indices */ strcpy( chLbl, chLineLbl(t) ); fprintf(ioPUN, "%li\t%s" , index , chLbl ); /* stat weights */ fprintf(ioPUN, "\t%.0f\t%.0f", (*t.Lo()).g() ,(*t.Hi()).g()); /* energy difference, gf */ fprintf(ioPUN, "\t%.2f\t%.3e", t.EnergyWN() ,t.Emis().gf()); fprintf(ioPUN, "\n"); } else { /* not first call, so do level populations and indices defined above */ if( (*t.Hi()).Pop() > xLimit ) { /* >>chng 05 may 08, add abundances, which for iso-seq species is * the density of the parent ion, for other lines, is unity. * had not been included so pops for iso seq were rel to parent ion. * caught by John Everett */ /* multiplication by abundance no longer necessary since iso pops now denormalized */ fprintf(ioPUN,"%li\t%.2e\t%.2e\n", index, (*t.Lo()).Pop(), (*t.Hi()).Pop() ); } } } /*SaveNewContinuum produce the 'save new continuum' output */ STATIC void SaveNewContinuum(FILE * ioPUN ) { long int ipLo, ipHi, j , ncells; double wllo, wlhi; double *energy, *cont_incid, *cont_atten, *diffuse_in, *diffuse_out, *emis_lines_out, *emis_lines_in; /* get the low limit */ wllo = rfield.anu[0]; /* get high-energy limit */ wlhi = rfield.anu[rfield.nflux-1]; /* use native continuum mesh */ ipLo = ipoint(wllo)-1; ipHi = ipoint(wlhi)-1; ncells = ipHi - ipLo + 1; /* now allocate the space */ energy = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ); cont_incid = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ); cont_atten = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ); diffuse_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ); diffuse_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ); emis_lines_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ); emis_lines_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ); /*emis_lines_pump_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double)); emis_lines_pump_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double));*/ /* now convert to rydbergs */ for( j=0; j 0. ) { SaveResults1Line(ioPUN,(char*)LineSv[i].chALab,LineSv[i].wavelength, LineSv[i].SumLine[0],"Line "); } } SaveResults1Line(ioPUN," ",0,0.,"Flush"); fprintf( ioPUN, " \n" ); fprintf( ioPUN, "BEGIN COLUMN DENSITIES\n" ); /* this dumps out the whole array,*/ /* following loop relies on LIMELM being 30, assert it here in case * this is ever changed */ ASSERT( LIMELM == 30 ); /* this order of indices is to keep 30 as the fastest variable, * and the 32 (LIMELM+1) as the slower one */ for( nelem=0; nelem>chng 02 apr 24, do not print type */ /* single column for input into data base */ if( strcmp(save.chPunRltType,"column") == 0 ) fprintf( ioPUN, "\n" ); } /* only put cr if we did not just put one */ if( strcmp(save.chPunRltType,"array ") == 0 ) fprintf( ioPUN, " \n" ); ipLine = 0; } } else if( strcmp(chFunction,"Flush") == 0 ) { if( ipLine > 0 ) { /* this is an option to print many emission lines across an output line, * the array option, or a single column of numbers, the column option * that appears on the "save results" and "save intensity" commands */ /* usual array 6 wide */ for( i=0; i < ipLine; i++ ) { fprintf( ioPUN, " %4.4s", chLabSave[i] ); prt_wl( ioPUN, wavelength[i] ); fprintf( ioPUN,"\t%.3e", xIntenSave[i]); /* >>chng 02 apr 24, do not print type */ /* single column for input into data base */ if( strcmp(save.chPunRltType,"column") == 0 ) fprintf( ioPUN, "\n" ); } if( strcmp(save.chPunRltType,"array ") == 0 ) fprintf( ioPUN, " \n" ); } } else { fprintf( ioQQQ, " SaveResults1Line called with insane request=%5.5s\n", chFunction ); cdEXIT(EXIT_FAILURE); } return; } static const int NENR_GAUNT = 37; static const int NTE_GAUNT = 21; /*SaveGaunts called by save gaunts command to output gaunt factors */ STATIC void SaveGaunts(FILE* ioPUN) { long int i, charge, ite, j; realnum ener[NENR_GAUNT], ste[NTE_GAUNT], z; double tempsave; double g[NENR_GAUNT][NTE_GAUNT], gfac; DEBUG_ENTRY( "SaveGaunts()" ); /* this routine is called from the PUNCH GAUNTS command * it drives the gaunt factor routine to save gaunts over full range */ tempsave = phycon.te; for( i=0; i < NTE_GAUNT; i++ ) { ste[i] = 0.5f*i; } for( i=0; i < NENR_GAUNT; i++ ) { ener[i] = 0.5f*i - 8.f; } for( charge = 1; charge 0 ) { sprintf( chStringHold, ", %f", grid.interpParameters[optimize.nOptimiz][j] ); strcat( chGridParam, chStringHold ); } fprintf( pnunit, "\t%f", grid.interpParameters[optimize.nOptimiz][j] ); } fprintf( pnunit, "\t%s\n", chGridParam ); }