00001
00002
00003
00004
00005
00006
00007 #include "cddefines.h"
00008 #include "dense.h"
00009 #include "deuterium.h"
00010 #include "struc.h"
00011 #include "rfield.h"
00012 #include "elementnames.h"
00013 #include "atoms.h"
00014 #include "physconst.h"
00015 #include "iterations.h"
00016 #include "pressure.h"
00017 #include "mole.h"
00018 #include "trace.h"
00019 #include "radius.h"
00020 #include "thermal.h"
00021 #include "heavy.h"
00022 #include "wind.h"
00023 #include "init.h"
00024 #include "iso.h"
00025 #include "h2.h"
00026 #include "co.h"
00027 #include "monitor_results.h"
00028 #include "taulines.h"
00029
00030
00031
00032
00033
00034 void InitSimPostparse( void )
00035 {
00036 DEBUG_ENTRY( "InitSimPostparse()" );
00037 static double one=1.0;
00038
00039 struc.nzonePreviousIteration = -1;
00040
00041 thermal.thist = 0.;
00042 thermal.tlowst = 1e20f;
00043 thermal.nUnstable = 0;
00044 thermal.lgUnstable = false;
00045
00046 colliders.init();
00047
00048 mole_global.init();
00049 mole_global.zero();
00050
00051 for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00052 {
00053 for( long ion=0; ion<nelem+2; ++ion )
00054 mole.set_location( nelem, ion, &(dense.xIonDense[nelem][ion]) );
00055 }
00056
00057 findspecieslocal("e-")->location = &(dense.eden);
00058 findspecieslocal("grn")->location = &(mole.grain_area);
00059 findspecieslocal("PHOTON")->location = &one;
00060 findspecieslocal("CRPHOT")->location = &one;
00061 findspecieslocal("CRP")->location = &one;
00062 if( deut.lgElmtOn )
00063 {
00064 findspecieslocal("D")->location = &(deut.xIonDense[0]);
00065 findspecieslocal("D+")->location = &(deut.xIonDense[1]);
00066 }
00067
00068
00069 mole_create_react();
00070
00071
00072
00073 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
00074 (*diatom)->H2_Reset();
00075
00076
00077 for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00078 for( long ion=0; ion<=nelem; ++ion )
00079 Heavy.xLyaHeavy[nelem][ion] = 0;
00080
00081
00082 if( radius.StopRadius[0] > 0. )
00083 {
00084 for( long j=0; j < iterations.iter_malloc; j++ )
00085 {
00086 if( radius.StopRadius[j] > radius.Radius )
00087 radius.StopThickness[j] = radius.StopRadius[j] - radius.Radius;
00088 else
00089 {
00090 fprintf( ioQQQ, " PROBLEM stopping radius is <= inner radius. Bailing out.\n" );
00091 cdEXIT(EXIT_FAILURE);
00092 }
00093 }
00094 }
00095
00096
00097 wind.DiskRadius = 0;
00098 if( wind.lgDisk )
00099 wind.DiskRadius = radius.Radius;
00100
00101 iterations.lgLastIt = false;
00102
00103 rfield_opac_zero( 0 , rfield.nupper );
00104
00105 rfield.lgUSphON = false;
00106
00107
00108 rfield.ipEnergyBremsThin = 0;
00109 rfield.EnergyBremsThin = 0.;
00110 rfield.comtot = 0.;
00111 rfield.cmcool = 0.;
00112 rfield.cinrat = 0.;
00113 rfield.extin_mag_B_point = 0.;
00114 rfield.extin_mag_V_point = 0.;
00115 rfield.extin_mag_B_extended = 0.;
00116 rfield.extin_mag_V_extended = 0.;
00117 rfield.EnerGammaRay = 7676.;
00118
00119 for( vector<TransitionList>::iterator it = AllTransitions.begin(); it != AllTransitions.end(); ++it )
00120 {
00121 for( TransitionList::iterator tr = it->begin(); tr != it->end(); ++tr )
00122 {
00123 tr->Emis().TauTrack().clear();
00124 }
00125 }
00126
00127
00128 wind.AccelGravity = (realnum)(GRAV_CONST*wind.comass*SOLAR_MASS/POW2(radius.Radius)*
00129 (1.-wind.DiskRadius/radius.Radius) );
00130 if( trace.lgTrace && trace.lgWind )
00131 {
00132 fprintf(ioQQQ," InitSimPostparse sets AccelGravity %.3e lgDisk?%c\n",
00133 wind.AccelGravity ,
00134 TorF(wind.lgDisk) );
00135 }
00136
00137
00138 pressure.PresInteg = 0.;
00139 pressure.PresIntegElecThin = 0.;
00140 pressure.pinzon = 0.;
00141
00142
00143 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00144 {
00145 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
00146 {
00147 iso_sp[ipISO][nelem].Reset();
00148
00149 if( nelem>ipHELIUM && !dense.lgElmtOn[nelem] )
00150 {
00151 iso_sp[ipISO][nelem].numLevels_max = 0;
00152 iso_sp[ipISO][nelem].nCollapsed_max = 0;
00153 iso_sp[ipISO][nelem].n_HighestResolved_max = 0;
00154
00155 iso_sp[ipISO][nelem].numLevels_local = 0;
00156 iso_sp[ipISO][nelem].nCollapsed_local = 0;
00157 iso_sp[ipISO][nelem].n_HighestResolved_local = 0;
00158 }
00159 else
00160 {
00161 iso_update_num_levels( ipISO, nelem );
00162
00163 ASSERT( iso_sp[ipISO][nelem].nCollapsed_max > 0 || iso_ctrl.lgCompileRecomb[ipISO] );
00164 }
00165 }
00166 }
00167
00168 if( iso_ctrl.lgPrintNumberOfLevels )
00169 {
00170 fprintf(ioQQQ,"\n\n Number of levels in ions treated by iso sequences.\n");
00171 fprintf(ioQQQ," ISO Element hi-n(l-resolved) #(l-resolved) n(collapsed)\n");
00172
00173 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00174 {
00175 for( long nelem=ipISO; nelem<LIMELM; ++nelem )
00176 {
00177
00178 fprintf(ioQQQ," %s %s %4li %4li %4li \n",
00179 iso_ctrl.chISO[ipISO] ,
00180 elementnames.chElementSym[nelem],
00181 iso_sp[ipISO][nelem].n_HighestResolved_max,
00182 iso_sp[ipISO][nelem].numLevels_max-iso_sp[ipISO][nelem].nCollapsed_max,
00183 iso_sp[ipISO][nelem].nCollapsed_max);
00184 }
00185 }
00186 }
00187 atoms.p2nit = 0.;
00188 atoms.d5200r = 0.;
00189
00190 MonitorResults.SumErrorCaseMonitor = 0.;
00191 MonitorResults.nSumErrorCaseMonitor = 0;
00192
00193 return;
00194 }