cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*RT_OTS compute diffuse fields due to H, He atoms, ion, triplets, metal recombination, 00004 * called by ConvBase */ 00005 /*RT_OTS_AddLine add local destruction of lines to ots field */ 00006 /*RT_OTS_AddCont add local destruction of continuum to ots field */ 00007 /*RT_OTS_Update sum flux, otscon, otslin, ConInterOut, outlin, to form SummeDif, SummedCon SummedOcc */ 00008 /*RT_OTS_Zero - zero out some vectors - 00009 * this is only called when code initialized by ContSetIntensity */ 00010 /*RT_OTS_ChkSum sanity check confirms summed continua reflect contents of individuals */ 00011 #include "cddefines.h" 00012 #include "taulines.h" 00013 #include "opacity.h" 00014 #include "dense.h" 00015 #include "iso.h" 00016 #include "hmi.h" 00017 #include "h2.h" 00018 #include "rfield.h" 00019 #include "conv.h" 00020 #include "rt.h" 00021 #include "mole_co_atom.h" 00022 #include "atomfeii.h" 00023 #include "heavy.h" 00024 #include "he.h" 00025 #include "trace.h" 00026 00027 /* this flag may be used for debugging ots rate changes */ 00028 static int nOTS_Line_type = -1; 00029 static int nOTS1=-1 , nOTS2=-1; 00030 /*add local destruction of continuum to ots field */ 00031 STATIC void RT_OTS_AddCont( 00032 /* the ots rate itself */ 00033 realnum ots, 00034 /* pointer to continuum cell for ots, on f scale */ 00035 long int ip); 00036 00037 /* =================================================================== */ 00038 00039 void RT_OTS(void) 00040 { 00041 long int 00042 ipla, 00043 ipISO , 00044 nelem, 00045 n; 00046 realnum 00047 difflya, 00048 esc, 00049 ots; 00050 00051 /* the Bowen HeII yield 00052 * >>chng 06 aug 08, from 0.3 to 0.4, update with netzer */ 00053 # define BOWEN 0.5f 00054 long int ipHi, 00055 ipLo; 00056 00057 double bwnfac; 00058 double ots660; 00059 realnum cont_phot_destroyed; 00060 double save_lya_dest, 00061 save_he2lya_dest; 00062 00063 double save_he2rec_dest; 00064 00065 /*static long int nCall=0; 00066 ++nCall; 00067 fprintf(ioQQQ,"debugggtos enter %li\n", nCall );*/ 00068 00069 DEBUG_ENTRY( "RT_OTS()" ); 00070 00071 /************************************************************************** 00072 * 00073 * the bowen HeII - OIII fluorescense problem 00074 * 00075 **************************************************************************/ 00076 nOTS_Line_type = 0; 00077 nelem = ipHELIUM; 00078 if( dense.lgElmtOn[nelem] ) 00079 { 00080 /* conversion per unit atom to OIII, at end of sub we divide by it, 00081 * to fix lines back to proper esc/dest probs */ 00082 bwnfac = BOWEN * MAX2(0.f,1.f- Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pesc - 00083 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pelec_esc ); 00084 00085 /* the last factor accounts for fact that two photons are produced, 00086 * and the branching ratio */ 00087 ots660 = Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Aul* 00088 StatesElem[ipH_LIKE][nelem][ipH2p].Pop*dense.xIonDense[nelem][nelem+1]* 00089 /*>>chng 06 aug 08, rm 0.8 factor from below, renorm aft discuss with Netzer */ 00090 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pdest *BOWEN*2.0; 00091 00092 /* now add this to the ots field */ 00093 if( ots660 > SMALLFLOAT ) 00094 RT_OTS_AddLine(ots660 , he.ip660 ); 00095 00096 /* decrease the destruction prob by the amount we will add elsewhere, 00097 * ok since dest probs updated on every iteration*/ 00098 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pdest *= (realnum)bwnfac; 00099 ASSERT( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pdest >= 0. ); 00100 { 00101 /* debugging code for line oscillation problems 00102 * most often Lya OTS oscillations*/ 00103 /*@-redef@*/ 00104 enum {DEBUG_LOC=false}; 00105 /*@+redef@*/ 00106 if( DEBUG_LOC ) 00107 { 00108 fprintf(ioQQQ,"DEBUG HeII Bowen nzone %li bwnfac:%.2e bwnfac esc:%.2e ots660 %.2e\n", 00109 nzone, 00110 bwnfac , 00111 bwnfac/BOWEN , 00112 ots660 ); 00113 } 00114 } 00115 } 00116 00117 else 00118 { 00119 bwnfac = 1.; 00120 } 00121 00122 /* save Lya loss rates so we can reset at end */ 00123 save_lya_dest = Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pdest; 00124 00125 /* >>chng 03 may 28, hydro.dstfe2lya should not go into the ots field since it was 00126 * actually lost in exciting FeII 00127 * this is not ready to be used - need to recover ots with small feii atom 00128 fprintf(ioQQQ,"feii debug\t%.3e\t%.3e\n", 00129 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest , hydro.dstfe2lya); 00130 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest = MAX2(0.f, 00131 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest-hydro.dstfe2lya);*/ 00132 00133 /* this is option to kill Lya ots rates, 00134 * rfield.lgLyaOTS is usually true (1), and set false (0) with 00135 * no lya ots command */ 00136 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pdest *= rfield.lgLyaOTS; 00137 00138 /* option to kill heii lya and rec continua ots */ 00139 save_he2lya_dest = Transitions[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Emis->Pdest; 00140 Transitions[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Emis->Pdest *= rfield.lgHeIIOTS; 00141 00142 /* option to kill heii lya and rec continua ots */ 00143 save_he2rec_dest = iso.RadRecomb[ipH_LIKE][ipHELIUM][ipH1s][ipRecRad]; 00144 iso.RadRecomb[ipH_LIKE][ipHELIUM][ipH1s][ipRecRad] *= rfield.lgHeIIOTS; 00145 00146 nOTS_Line_type = 1; 00147 00148 /* make ots fields due to lines and continua of species treated with unified 00149 * isoelectronic sequence */ 00150 /* loop over all elements */ 00151 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00152 { 00153 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00154 { 00155 nOTS1 = ipISO; 00156 nOTS2 = nelem; 00157 /* do if this stage exists */ 00161 if( (dense.IonHigh[nelem] >= nelem+1-ipISO ) ) 00162 { 00163 /* generate line ots rates */ 00164 /* now loop over all possible levels, but cannot include two photon 00165 * since there is no pointer to this continuum */ 00166 /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */ 00167 for( ipHi=1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ ) 00168 { 00169 for( ipLo=0; ipLo < ipHi; ipLo++ ) 00170 { 00171 /* this signifies a fake line */ 00172 /* >>chng 03 may 19, DEST0 is the smallest possible 00173 * dest prob - not a real number, don't add to ots field */ 00174 if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont < 1 || 00175 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pdest<= DEST0 ) 00176 continue; 00177 00178 /* ots rates, the destp prob was set in hydropesc */ 00179 Transitions[ipISO][nelem][ipHi][ipLo].Emis->ots = 00180 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul* 00181 StatesElem[ipISO][nelem][ipHi].Pop*dense.xIonDense[nelem][nelem+1-ipISO]* 00182 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pdest; 00183 00184 ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Emis->ots >= 0. ); 00185 /* way to kill lyalpha ots rates 00186 if( nelem==ipHYDROGEN && ipHi==ipH2p && ipLo==ipH1s ) 00187 { 00188 Transitions[ipISO][nelem][ipHi][ipLo].Emis->ots = 0.; 00189 } */ 00190 00191 /* finally dump the ots rate into the stack */ 00192 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->ots > SMALLFLOAT ) 00193 RT_OTS_AddLine(Transitions[ipISO][nelem][ipHi][ipLo].Emis->ots, 00194 Transitions[ipISO][nelem][ipHi][ipLo].ipCont ); 00195 } 00196 } 00197 { 00198 /* debugging code for line oscillation problems 00199 * most often Lya OTS oscillations*/ 00200 /*@-redef@*/ 00201 enum {DEBUG_LOC=false}; 00202 /*@+redef@*/ 00203 if( DEBUG_LOC ) 00204 { 00205 long ip; 00206 if( ipISO==0 && nelem==0 && nzone>500 ) 00207 { 00208 ipHi = 2; 00209 ipLo = 0; 00210 ip = Transitions[ipISO][nelem][ipHi][ipLo].ipCont; 00211 fprintf(ioQQQ,"DEBUG hlyaots\t%.2f\tEdenTrue\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 00212 fnzone, 00213 dense.EdenTrue , 00214 Transitions[ipISO][nelem][ipHi][ipLo].Emis->ots, 00215 opac.opacity_abs[ip-1], 00216 StatesElem[ipISO][nelem][ipHi].Pop*dense.xIonDense[nelem][nelem+1-ipISO], 00217 StatesElem[ipISO][nelem][ipHi].Pop, 00218 dense.xIonDense[nelem][nelem+1-ipISO], 00219 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pdest, 00220 rfield.otslinNew[ip-1], 00221 rfield.otslin[ip-1]); 00222 } 00223 } 00224 } 00225 00226 /************************************************************************** 00227 * 00228 * ots recombination bound-free b-f continua continuum 00229 * 00230 **************************************************************************/ 00231 00232 /* put in OTS continuum */ 00233 /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */ 00234 for( n=0; n < iso.numLevels_local[ipISO][nelem]; n++ ) 00235 { 00236 cont_phot_destroyed = (realnum)(iso.RadRecomb[ipISO][nelem][n][ipRecRad]* 00237 (1. - iso.RadRecomb[ipISO][nelem][n][ipRecEsc])*dense.eden* 00238 dense.xIonDense[nelem][nelem+1-ipISO]); 00239 ASSERT( cont_phot_destroyed >= 0. ); 00240 00241 /* continuum energy index used in this routine is decremented by one there */ 00242 RT_OTS_AddCont(cont_phot_destroyed,iso.ipIsoLevNIonCon[ipISO][nelem][n]); 00243 /* debugging code for rec continua */ 00244 { 00245 /*@-redef@*/ 00246 enum {DEBUG_LOC=false}; 00247 /*@+redef@*/ 00248 if( DEBUG_LOC && nzone > 400 && nelem==0 && n==2 ) 00249 { 00250 long ip = iso.ipIsoLevNIonCon[ipISO][nelem][n]-1; 00251 fprintf(ioQQQ,"hotsdebugg\t%.3f\t%li\th con ots\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 00252 fnzone, 00253 n , 00254 StatesElem[ipISO][ipHYDROGEN][2].Pop*dense.xIonDense[ipHYDROGEN][1], 00255 cont_phot_destroyed, 00256 cont_phot_destroyed/opac.opacity_abs[ip], 00257 rfield.otsconNew[ip] , 00258 rfield.otscon[ip] , 00259 opac.opacity_abs[ip] , 00260 hmi.Hmolec[ipMHm] , 00261 hmi.HMinus_photo_rate); 00262 } 00263 } 00264 } 00265 } 00266 } 00267 } 00268 /* more debugging code for rec continua */ 00269 { 00270 /*@-redef@*/ 00271 enum {DEBUG_LOC=false}; 00272 /*@+redef@*/ 00273 if( DEBUG_LOC ) 00274 { 00275 nelem = 0; 00276 fprintf(ioQQQ,"hotsdebugg %li \t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 00277 nzone, 00278 rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][0]-1], 00279 rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][1]-1], 00280 rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][3]-1], 00281 rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][4]-1], 00282 rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][5]-1], 00283 rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][6]-1], 00284 opac.opacity_abs[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][6]-1]); 00285 } 00286 } 00287 00288 /* now reset Lya dest prob in case is was clobbered */ 00289 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pdest = (realnum)save_lya_dest; 00290 Transitions[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Emis->Pdest = (realnum)save_he2lya_dest; 00291 iso.RadRecomb[ipH_LIKE][ipHELIUM][ipH1s][ipRecRad] = save_he2rec_dest; 00292 00293 nelem = ipHELIUM; 00294 if( dense.lgElmtOn[nelem] && bwnfac > 0. ) 00295 { 00296 /* increase the destruction prob by the amount we decreased it above */ 00297 Transitions[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Emis->Pdest /= (realnum)bwnfac; 00298 } 00299 00300 if( trace.lgTrace ) 00301 { 00302 fprintf(ioQQQ," RT_OTS Pdest %.2e ots rate %.2e in otslinNew %.2e con opac %.2e\n", 00303 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pdest , Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->ots , 00304 rfield.otslinNew[Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1] , 00305 opac.opacity_abs[Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1] 00306 ); 00307 } 00308 00309 nOTS_Line_type = 2; 00310 /* recombination Lya for all elements not yet converted into std isoelectronc form */ 00311 for( nelem=NISO; nelem < LIMELM; nelem++ ) 00312 { 00313 long int ion; 00314 /* do not include species treated in iso-electronic fashon in the following, 00315 * these were treated above */ 00316 for( ion=0; ion < nelem+1-NISO; ion++ ) 00317 { 00318 if( dense.xIonDense[nelem][ion+1] > 0. ) 00319 { 00320 nOTS1 = nelem; 00321 nOTS2 = ion; 00322 /* now do the recombination Lya */ 00323 ipla = Heavy.ipLyHeavy[nelem][ion]; 00324 ASSERT( ipla>0 ); 00325 esc = opac.ExpmTau[ipla-1]; 00326 /* xLyaHeavy is set to a fraction of total rad rec in ion_recomb, includes eden */ 00327 difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1]; 00328 /* >>chng 00 dec 22, from MIN2 to MAX2, MIN2 had effect of always 00329 * setting the ots rates to zero */ 00330 ots = difflya*MAX2(0.f,1.f-esc); 00331 /*if( nelem==6 && ion==2 ) 00332 fprintf(ioQQQ," debugggnly\t %.2f\t%.2e\n",fnzone, ots );*/ 00333 ASSERT( ots >= 0.); 00334 /*if( iteration == 2 && nzone>290 && ipla== 2339 ) 00335 fprintf(ioQQQ,"recdebugg1 %.2e %li %li %.2e %.2e \n", 00336 ots, nelem, ion, 00337 esc , dense.xIonDense[nelem][ion+1]);*/ 00338 if( ots > SMALLFLOAT ) 00339 RT_OTS_AddLine(ots,ipla); 00340 00341 /* now do the recombination balmer lines */ 00342 ipla = Heavy.ipBalHeavy[nelem][ion]; 00343 esc = opac.ExpmTau[ipla-1]; 00344 /* xLyaHeavy is set to a fraction of total rad rec in ion_recomb, includes eden */ 00345 difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1]; 00346 /* >>chng 00 dec 22, from MIN2 to MAX2, MIN2 had effect of always 00347 * setting the ots rates to zero */ 00348 ots = difflya*MAX2(0.f,1.f-esc); 00349 ASSERT( ots >= 0.); 00350 /*if( iteration == 2 &&nzone==294 && ipla== 2339 ) 00351 fprintf(ioQQQ,"recdebugg2 %.2e %li %li\n", ots, nelem, ion );*/ 00352 if( ots > SMALLFLOAT ) 00353 RT_OTS_AddLine(ots,ipla); 00354 } 00355 } 00356 } 00357 00358 nOTS_Line_type = 3; 00359 /* do OTS and outward parts of FeII lines, if large atom is turned on */ 00360 FeII_OTS(); 00361 00362 nOTS_Line_type = 4; 00363 /* do the main set of level1 lines */ 00364 # if 0 00365 { 00366 # include "lines_service.h" 00367 if( nzone>290 ) DumpLine( &TauLines[74] ); 00368 } 00369 # endif 00370 /*# define FRACNEW 1.0f 00371 # undef FRACNEW*/ 00372 for( nOTS1=1; nOTS1 < nLevel1; nOTS1++ ) 00373 { 00374 /*realnum otsold = TauLines[nOTS1].ots;*/ 00375 TauLines[nOTS1].Emis->ots = TauLines[nOTS1].Hi->Pop * TauLines[nOTS1].Emis->Aul * TauLines[nOTS1].Emis->Pdest; 00376 /* * TauLines[nOTS1].ColOvTot;*/ 00377 /*if( nzone ) 00378 TauLines[nOTS1].ots = TauLines[nOTS1].ots*FRACNEW + otsold*(1.f-FRACNEW);*/ 00379 if( TauLines[nOTS1].Emis->ots > SMALLFLOAT ) 00380 RT_OTS_AddLine( TauLines[nOTS1].Emis->ots , TauLines[nOTS1].ipCont); 00381 } 00382 00383 nOTS_Line_type = 5; 00384 /* do the level2 level 2 lines */ 00385 for( nOTS1=0; nOTS1 < nWindLine; nOTS1++ ) 00386 { 00387 if( TauLine2[nOTS1].Hi->IonStg < TauLine2[nOTS1].Hi->nelem+1-NISO ) 00388 { 00389 TauLine2[nOTS1].Emis->ots = TauLine2[nOTS1].Hi->Pop * TauLine2[nOTS1].Emis->Aul * TauLine2[nOTS1].Emis->Pdest; 00390 if( TauLine2[nOTS1].Emis->ots > SMALLFLOAT ) 00391 RT_OTS_AddLine( TauLine2[nOTS1].Emis->ots , TauLine2[nOTS1].ipCont); 00392 } 00393 } 00394 /*fprintf(ioQQQ,"debuggne1 \t%.2e\t%.2e\t%.2e\t%.2e\n", 00395 TauLine2[743].ots/SDIV(opac.opacity_abs[743]), 00396 TauLine2[744].ots/SDIV(opac.opacity_abs[744]) , 00397 TauLine2[743].ots,opac.opacity_abs[743]);*/ 00398 /*CO_OTS - add CO lines to ots fields, called by RT_OTS */ 00399 00400 nOTS_Line_type = 6; 00401 CO_OTS(); 00402 00403 /* the large H2 molecule */ 00404 nOTS_Line_type = 7; 00405 H2_RT_OTS(); 00406 /*fprintf(ioQQQ,"DEBUG ggtos exit\n");*/ 00407 return; 00408 } 00409 00410 /* =================================================================== */ 00411 00412 void RT_OTS_AddLine(double ots, 00413 /* pointer on the f scale */ 00414 long int ip ) 00415 { 00416 00417 DEBUG_ENTRY( "RT_OTS_AddLine()" ); 00418 00419 /* add ots due to line destruction to radiation field */ 00420 00421 /* return if outside bounds of this continuum source, ip > rfield.nflux 00422 * first case ip==0 happens when called with dummy line */ 00423 if( ip==0 || ip > rfield.nflux ) 00424 { 00425 return; 00426 } 00427 00428 /*the local ots rate must be non-negative */ 00429 ASSERT( ots >= 0. ); 00430 /* continuum pointer must be positive */ 00431 ASSERT( ip > 0 ); 00432 00433 /* add locally destroyed flux of photons to line OTS array */ 00434 /* check whether local gas opacity (units cm-1) is positive, if so 00435 * convert line destruction rate into ots rate by dividing by it */ 00436 if( opac.opacity_abs[ip-1] > 0. ) 00437 { 00438 rfield.otslinNew[ip-1] += (realnum)(ots/opac.opacity_abs[ip-1]); 00439 } 00440 /* first iteration is 1, second is two */ 00441 { 00442 /*@-redef@*/ 00443 enum {DEBUG_LOC=false}; 00444 /*@+redef@*/ 00445 if( DEBUG_LOC && /*iteration == 2 && nzone>294 &&*/ ip== 2363 /*&& ots > 1e16*/) 00446 { 00447 fprintf(ioQQQ,"DEBUG ots, opc, otsr %.3e\t%.3e\t%.3e\t", 00448 ots , 00449 opac.opacity_abs[ip-1], 00450 ots/opac.opacity_abs[ip-1] ); 00451 fprintf(ioQQQ,"iteration %li type %i %i %i \n", 00452 iteration, 00453 nOTS_Line_type, 00454 nOTS1,nOTS2 ); 00455 } 00456 } 00457 return; 00458 } 00459 00460 /* =================================================================== */ 00461 00462 /*add local destruction of continuum to ots field */ 00463 STATIC void RT_OTS_AddCont( 00464 /* the ots rate itself */ 00465 realnum ots, 00466 /* pointer to continuum cell for ots, on f scale */ 00467 long int ip) 00468 { 00469 00470 DEBUG_ENTRY( "RT_OTS_AddCont()" ); 00471 00472 /* 00473 * routine called to add ots due to continuum destruction to 00474 * radiation field 00475 */ 00476 00477 /* check if outside bounds of this continuum source */ 00478 if( ip > rfield.nflux ) 00479 { 00480 return; 00481 } 00482 00483 ASSERT( ip> 0 ); 00484 ASSERT( ots >= 0. ); 00485 ASSERT( ip <= rfield.nupper ); 00486 00487 /* add locally destroyed flux of photons to continuum OTS array */ 00488 /* check whether local gas opacity (units cm-1) is positive, if so 00489 * convert continuum destruction rate into ots rate by dividing by it */ 00490 if( opac.opacity_abs[ip-1] > 0. ) 00491 { 00492 rfield.otsconNew[ip-1] += (realnum)(ots/opac.opacity_abs[ip-1]); 00493 /*if( ots>0. ) fprintf(ioQQQ, 00494 "buggg ip %li ots %.2e opac %.2e cont_phot_destroyed %.2e\n", 00495 ip, ots , opac.opacity_abs[ip-1] ,rfield.otsconNew[ip-1]);*/ 00496 } 00497 return; 00498 } 00499 00500 #if 0 00501 /* check whether Lya OTS rate is oscillating */ 00502 STATIC bool lgLyaOTS_oscil( realnum ots_new ) 00503 { 00504 static realnum old , older; 00505 bool lgReturn = false; 00506 00507 if( !conv.nTotalIoniz || conv.lgSearch ) 00508 { 00509 /* this is very first call on this iteration, zero everything out */ 00510 old = 0.; 00511 older = 0.; 00512 } 00513 00514 /* is sign of change oscillating? */ 00515 if( conv.nPres2Ioniz>1 && (ots_new-old) * (old-older) < 0.) 00516 lgReturn = true; 00517 00518 older = old; 00519 old = ots_new; 00520 return( lgReturn ); 00521 } 00522 00523 /* >>chng 04 jan 26, check whether Lya is oscillating */ 00524 if( lgLyaOTS_oscil( rfield.otslinNew[Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1] ) ) 00525 { 00526 fprintf(ioQQQ,"DEBUG lya ots %.2f %.3e\n", 00527 fnzone, 00528 rfield.otslinNew[Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1] ); 00529 } 00530 #endif 00531 /* =================================================================== */ 00532 00533 /*RT_OTS_Update update ots rates, called in ConvBase */ 00534 void RT_OTS_Update( 00535 /* summed ots rates */ 00536 double *SumOTS , 00537 /* array index for constituent that changed the most 00538 long int *ipOTSchange, */ 00539 /* limit on how large a change to allow, 0 for no limit */ 00540 double BigFrac) 00541 { 00542 long int i; 00543 00544 static bool lgNeedSpace=true; 00545 static double *old_OTS_line_x_opac , *old_OTS_cont_x_opac; 00546 realnum FacBig , FacSml; 00547 00548 DEBUG_ENTRY( "RT_OTS_Update()" ); 00549 00550 if( BigFrac <= 0. ) 00551 BigFrac = 10.; 00552 FacBig = (realnum)(1. + BigFrac); 00553 FacSml = 1.f/FacBig; 00554 00555 /* create space to save arrays if needed, one time per coreload */ 00556 if( lgNeedSpace ) 00557 { 00558 old_OTS_line_x_opac = (double*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(double)) ); 00559 old_OTS_cont_x_opac = (double*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(double)) ); 00560 } 00561 lgNeedSpace = false; 00562 00563 *SumOTS = 0.; 00564 00565 /* option to kill ots rates with no ots lines command */ 00566 if( rfield.lgKillOTSLine ) 00567 { 00568 for( i=0; i < rfield.nflux; i++ ) 00569 { 00570 rfield.otslinNew[i] = 0.; 00571 } 00572 } 00573 00574 /* during search phase set ots to current values */ 00575 if( nzone==0 ) 00576 { 00577 for( i=0; i < rfield.nflux; i++ ) 00578 { 00579 old_OTS_line_x_opac[i] = rfield.otslinNew[i] * opac.opacity_abs[i]; 00580 old_OTS_cont_x_opac[i] = rfield.otsconNew[i] * opac.opacity_abs[i]; 00581 } 00582 } 00583 00584 /* remember largest change in ots rates */ 00585 *SumOTS = 0.; 00586 /* now update new ots rates */ 00587 for( i=0; i < rfield.nflux; ++i ) 00588 { 00589 double OTS_line_x_opac = rfield.otslinNew[i] * opac.opacity_abs[i]; 00590 double OTS_cont_x_opac = rfield.otsconNew[i] * opac.opacity_abs[i]; 00591 double CurrentInverseOpacity = 1./MAX2( SMALLDOUBLE , opac.opacity_abs[i] ); 00592 00593 /* get new ots line rates for each cell, but don't let change be too big */ 00594 if( OTS_line_x_opac > old_OTS_line_x_opac[i]*FacBig ) 00595 { 00596 rfield.otslin[i] = rfield.otslin[i]*FacBig; 00597 } 00598 else if( OTS_line_x_opac < old_OTS_line_x_opac[i]*FacSml ) 00599 { 00600 rfield.otslin[i] = rfield.otslin[i]*FacSml; 00601 } 00602 else 00603 { 00604 rfield.otslin[i] = rfield.otslinNew[i]; 00605 } 00606 00607 /* get new ots continuum rates for each cell */ 00608 if( OTS_cont_x_opac > old_OTS_cont_x_opac[i]*FacBig ) 00609 { 00610 rfield.otscon[i] = rfield.otscon[i]*FacBig; 00611 } 00612 else if( OTS_cont_x_opac < old_OTS_cont_x_opac[i]*FacSml ) 00613 { 00614 rfield.otscon[i] = rfield.otscon[i]*FacSml; 00615 } 00616 else 00617 { 00618 rfield.otscon[i] = rfield.otsconNew[i]; 00619 } 00620 00621 /* zero out the new otscon vector*/ 00622 rfield.otsconNew[i] = 0.; 00623 00624 /* zero out the new otslin vector*/ 00625 rfield.otslinNew[i] = 0.; 00626 00627 /* now update to new values */ 00628 old_OTS_line_x_opac[i] = rfield.otslin[i] * opac.opacity_abs[i]; 00629 old_OTS_cont_x_opac[i] = rfield.otscon[i] * opac.opacity_abs[i]; 00630 00631 /* this is local ots continuum created by destroyed diffuse continua, 00632 * currently only two-photon */ 00633 rfield.ConOTS_local_OTS_rate[i] = (realnum)((double)rfield.ConOTS_local_photons[i]*CurrentInverseOpacity); 00634 00635 /*if( i==14 ) 00636 { 00637 fprintf(ioQQQ,"DEBUG OTS rate %.2e %.2e %.2e %.2e \n", *SumOTS , 00638 rfield.otscon[14] , rfield.otslin[14],opac.opacity_abs[14]); 00639 fflush(ioQQQ); 00640 }*/ 00641 00642 /* remember sum of ots rates for convergence criteria */ 00643 *SumOTS += (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i]; 00644 00645 rfield.SummedDif[i] = rfield.otscon[i] + rfield.otslin[i] + rfield.outlin_noplot[i]+ 00646 rfield.ConInterOut[i]*rfield.lgOutOnly + rfield.outlin[i] + 00647 rfield.ConOTS_local_OTS_rate[i]; 00648 00649 rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i]; 00650 rfield.SummedOcc[i] = rfield.SummedCon[i]*rfield.convoc[i]; 00651 00652 } 00653 00654 00655 /* >>chng 02 oct 18, add this */ 00656 /* sum of accumulated flux from particular frequency to infinity */ 00657 rfield.flux_accum[rfield.nflux-1] = 0; 00658 for( i=1; i < rfield.nflux; i++ ) 00659 { 00660 rfield.flux_accum[rfield.nflux-i-1] = rfield.flux_accum[rfield.nflux-i] + 00661 rfield.SummedCon[rfield.nflux-i-1]; 00662 } 00663 /* this has to be positive since is sum of all photons in SummedCon */ 00664 ASSERT( rfield.flux_accum[0]> 0. ); 00665 00666 /* >>chng 02 jul 23, set to black body at local temp if in optically thick continuum, 00667 * between plasma frequency and energy where brems is thin */ 00668 ASSERT(rfield.ipPlasma>0 ); 00669 00670 /* >>chng 02 jul 25, set all radiation fields to zero below plasma frequency */ 00671 for( i=0; i < rfield.ipPlasma-1; i++ ) 00672 { 00673 rfield.otscon[i] = 0.; 00674 rfield.ConOTS_local_OTS_rate[i] = 0.; 00675 rfield.ConOTS_local_photons[i] = 0.; 00676 rfield.otslin[i] = 0.; 00677 rfield.SummedDif[i] = 0.; 00678 rfield.OccNumbBremsCont[i] = 0.; 00679 rfield.SummedCon[i] = 0.; 00680 rfield.SummedOcc[i] = 0.; 00681 rfield.ConInterOut[i] = 0.; 00682 } 00683 /* this loop evaluates occupation number for brems continuum, 00684 * only used for induced two photon emission */ 00685 if( rfield.ipEnergyBremsThin > 0 ) 00686 { 00687 for( i=rfield.ipPlasma-1; i < rfield.nflux; i++ ) 00688 { 00689 /* this corrects for opacity / optical depth in brems - brems opacity goes as 00690 * energy squared. */ 00691 /* when totally optically thin to brems rfield.ipEnergyBremsThin is zero, 00692 * so need this max */ 00693 realnum factor = MIN2(1.f,rfield.anu2[MAX2(0,rfield.ipEnergyBremsThin-1)] / rfield.anu2[i]); 00694 00695 /* occupation number for black body is 1/ (exp hn/kT) -1) */ 00696 rfield.OccNumbBremsCont[i] = (realnum)(1./(1./SDIV(rfield.ContBoltz[i]) - 1.)) * factor; 00697 } 00698 } 00699 return; 00700 00701 # if 0 00702 OTSOld = 0.; 00703 OTSNew = 0.; 00704 BigOTSNew = 0.; 00705 /* find summed ots rates, old and new */ 00706 for( i=0; i < rfield.nflux; i++ ) 00707 { 00708 OTSOld += rfield.otscon[i] + rfield.otslin[i]; 00709 OTSNew += rfield.otsconNew[i] + rfield.otslinNew[i]; 00710 if( BigOTSNew < rfield.otsconNew[i] + rfield.otslinNew[i] ) 00711 { 00712 BigOTSNew = rfield.otsconNew[i] + rfield.otslinNew[i]; 00713 } 00714 } 00715 00716 /* we now have old and new rates, what is the ratio, and by how much will 00717 * we allow this to change? */ 00718 if( BigFrac == 0. || conv.lgSearch || OTSOld < SMALLFLOAT ) 00719 { 00720 /* this branch, want a clear update of ots rates, either because 00721 * requested with BigFrac == 0, or we are in search phase */ 00722 FracOld = 0.; 00723 } 00724 else 00725 { 00726 if( OTSNew > OTSOld ) 00727 { 00728 chng = fabs(1. - OTSOld/SDIV(OTSNew) ); 00729 } 00730 else 00731 { 00732 chng = fabs(1. - OTSNew/SDIV(OTSOld) ); 00733 } 00734 00735 if( chng < BigFrac ) 00736 { 00737 /* this branch, old and new do not differ by much, so use new */ 00738 FracOld = 0.; 00739 } 00740 else 00741 { 00742 /* this branch, too large a difference between old and new, cap it to BigFrac */ 00743 FracOld = (1. - BigFrac / chng); 00744 ASSERT( FracOld >= 0. ); 00745 FracOld = MIN2( 0.25 , FracOld ); 00746 } 00747 } 00748 00749 /* fraction old and new ots rates */ 00750 FracNew = 1. - FracOld; 00751 fprintf(ioQQQ," DEBUG FracNew\t%.2e\n", FracNew); 00752 00753 /* remember largest change in ots rates */ 00754 changeOTS = 0.; 00755 /*fprintf(ioQQQ," sumcontinuum zone%li 1168=%e\n", nzone,rfield.otslin[1167]);*/ 00756 00757 for( i=0; i < rfield.nflux; i++ ) 00758 { 00759 /* >>chng 01 feb 01, define inverse opacity in safe manner */ 00760 double CurrentInverseOpacity = 1./MAX2( SMALLDOUBLE , opac.opacity_abs[i] ); 00761 00762 /* remember which energy had largest change in ots rates */ 00763 if( fabs( rfield.otscon[i]+rfield.otslin[i]-rfield.otsconNew[i]-rfield.otslinNew[i])> changeOTS) 00764 { 00765 changeOTS = 00766 fabs( rfield.otscon[i]+rfield.otslin[i]-rfield.otsconNew[i]-rfield.otslinNew[i]); 00767 } 00768 00769 /* >>chng 01 apr 10, only change ots with means if FracOld not zero, 00770 * this is to avoid taking means when this routine called by StartEndIter 00771 * to reset sums of rates */ 00772 if( BigFrac > 0. && conv.nTotalIoniz > 0 ) 00773 { 00774 /* the New vectors are the ones updated by the AddOTS routines, 00775 * and have the current OTS rates. The otscon and otslin vectors 00776 * have the rates from the previous refresh of the New vectors*/ 00777 /* here is the refresh. all were initially set to zero in call to 00778 * RT_OTS_Zero (below) from zero */ 00779 rfield.otscon[i] = (realnum)((rfield.otscon[i]*FracOld +rfield.otsconNew[i]*FracNew)); 00780 00781 /* this is local ots continuum created by destroyed diffuse continua, 00782 * currently only two-photon */ 00783 /* >>chng 00 dec 27, define this continuum and do not count 2-photon in old var */ 00784 rfield.ConOTS_local_OTS_rate[i] = (realnum)(rfield.ConOTS_local_photons[i]*CurrentInverseOpacity); 00785 00786 /* current ots will be average of old and new */ 00787 rfield.otslin[i] = (realnum)((rfield.otslin[i]*FracOld + rfield.otslinNew[i]*FracNew)); 00788 } 00789 00790 /* zero out the new otscon vector*/ 00791 rfield.otsconNew[i] = 0.; 00792 00793 /* zero out the new otslin vector*/ 00794 rfield.otslinNew[i] = 0.; 00795 00796 /* remember sum of ots rates for convergence criteria */ 00797 *SumOTS += (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i]; 00798 { 00799 /* following should be set true to print strongest ots contributors */ 00800 /*@-redef@*/ 00801 enum {DEBUG_LOC=false}; 00802 /*@+redef@*/ 00803 if( DEBUG_LOC && (nzone==378) ) 00804 { 00805 if( conv.nPres2Ioniz > 3 ) 00806 cdEXIT( EXIT_FAILURE ); 00807 fprintf(ioQQQ,"rtotsbugggg\t%li\t%.3e\t%.3e\t%.3e\t%.3e\n", 00808 conv.nPres2Ioniz, 00809 rfield.anu[i], 00810 opac.opacity_abs[i], 00811 rfield.otscon[i], 00812 rfield.otslin[i]); 00813 } 00814 } 00815 00816 /* >>chng 00 dec 27, include ConOTS_local_OTS_rate */ 00817 /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */ 00818 rfield.SummedDif[i] = rfield.otscon[i] + rfield.otslin[i] + rfield.outlin_noplot[i]+ 00819 rfield.ConInterOut[i]*rfield.lgOutOnly + rfield.outlin[i] + 00820 rfield.ConOTS_local_OTS_rate[i]; 00821 00822 rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i]; 00823 rfield.SummedOcc[i] = rfield.SummedCon[i]*rfield.convoc[i]; 00824 } 00825 00826 ASSERT(*SumOTS >= 0. ); 00827 00828 /* >>chng 02 oct 18, add this */ 00829 /* sum of accumulated flux from particular frequency to infinity */ 00830 rfield.flux_accum[rfield.nflux-1] = 0; 00831 for( i=1; i < rfield.nflux; i++ ) 00832 { 00833 rfield.flux_accum[rfield.nflux-i-1] = rfield.flux_accum[rfield.nflux-i] + 00834 rfield.SummedCon[rfield.nflux-i-1]; 00835 } 00836 /* this has to be positive since is sum of all photons in SummedCon */ 00837 ASSERT( rfield.flux_accum[0]> 0. ); 00838 00839 /* >>chng 02 jul 23, set to black body at local temp if in optically thick continuum, 00840 * between plasma frequency and energy where brems is thin */ 00841 ASSERT(rfield.ipPlasma>0 ); 00842 00843 /* >>chng 02 jul 25, set all radiation fields to zero below plasma frequency */ 00844 for( i=0; i < rfield.ipPlasma-1; i++ ) 00845 { 00846 rfield.otscon[i] = 0.; 00847 rfield.ConOTS_local_OTS_rate[i] = 0.; 00848 rfield.ConOTS_local_photons[i] = 0.; 00849 rfield.otslin[i] = 0.; 00850 rfield.SummedDif[i] = 0.; 00851 rfield.OccNumbBremsCont[i] = 0.; 00852 rfield.SummedCon[i] = 0.; 00853 rfield.SummedOcc[i] = 0.; 00854 rfield.ConInterOut[i] = 0.; 00855 } 00856 /* this loop evaluates occupation number for brems continuum, 00857 * only used for induced two photon emission */ 00858 if( rfield.ipEnergyBremsThin > 0 ) 00859 { 00860 for( i=rfield.ipPlasma-1; i < rfield.nflux; i++ ) 00861 { 00862 /* this corrects for opacity / optical depth in brems - brems opacity goes as 00863 * energy squared. */ 00864 /* when totally optically thin to brems rfield.ipEnergyBremsThin is zero, 00865 * so need this max */ 00866 realnum factor = MIN2(1.f,rfield.anu2[MAX2(0,rfield.ipEnergyBremsThin-1)] / rfield.anu2[i]); 00867 00868 /* occupation number for black body is 1/ (exp hn/kT) -1) */ 00869 rfield.OccNumbBremsCont[i] = (realnum)(1./(1./SDIV(rfield.ContBoltz[i]) - 1.)) * factor; 00870 } 00871 } 00872 00873 { 00874 /* following should be set true to print strongest ots contributors */ 00875 /*@-redef@*/ 00876 enum {DEBUG_LOC=false}; 00877 /*@+redef@*/ 00878 if( DEBUG_LOC && (nzone>0) ) 00879 { 00880 double BigOTS; 00881 int ipOTS=0; 00882 BigOTS = -1.; 00883 /* find the biggest ots contributor */ 00884 /*for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]; i < rfield.nflux; i++ )*/ 00885 for( i=0; i < rfield.nflux; i++ ) 00886 { 00887 if( (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i] > BigOTS ) 00888 { 00889 BigOTS = (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i]; 00890 ipOTS = (int)i; 00891 } 00892 } 00893 fprintf(ioQQQ, 00894 " sumcontinuunots zone %li SumOTS %.2e %.2eRyd opc:%.2e OTS%.2e lin%.2e con:%.2e %s %s cell%i FracNew %.2f \n", 00895 nzone, 00896 *SumOTS, 00897 rfield.anu[ipOTS] , 00898 opac.opacity_abs[ipOTS], 00899 BigOTS , 00900 rfield.otslin[ipOTS], 00901 rfield.otscon[ipOTS] , 00902 /* line label */ 00903 rfield.chLineLabel[ipOTS] , 00904 /* cont label*/ 00905 rfield.chContLabel[ipOTS] , 00906 ipOTS, 00907 FracNew); 00908 } 00909 } 00910 return; 00911 # endif 00912 } 00913 00914 /* =================================================================== */ 00915 00916 /*RT_OTS_Zero zero out some vectors - this is only called when code 00917 * initialized by ContSetIntensity */ 00918 void RT_OTS_Zero( void ) 00919 { 00920 long int i; 00921 00922 DEBUG_ENTRY( "RT_OTS_Zero()" ); 00923 00924 /* this loop goes up to nflux itself (<=) since the highest cell 00925 * will be used to pass unity through the code to verify integrations */ 00926 for( i=0; i <= rfield.nflux; i++ ) 00927 { 00928 rfield.SummedDif[i] = 0.; 00929 /* the four main ots vectors */ 00930 rfield.otscon[i] = 0.; 00931 rfield.otslin[i] = 0.; 00932 rfield.otslinNew[i] = 0.; 00933 rfield.otsconNew[i] = 0.; 00934 /* >>chng 05 mar 03, add this */ 00935 rfield.trans_coef_zone[i] = 1.f; 00936 rfield.trans_coef_total[i] = 1.f; 00937 00938 rfield.ConInterOut[i] = 0.; 00939 rfield.outlin[i] = 0.; 00940 rfield.outlin_noplot[i] = 0.; 00941 rfield.SummedDif[i] = 0.; 00942 /* "zero" for the summed con will be just the incident radiation field */ 00943 rfield.SummedCon[i] = rfield.flux[i]; 00944 rfield.SummedOcc[i] = rfield.SummedCon[i]*rfield.convoc[i]; 00945 rfield.ConOTS_local_photons[i] = 0.; 00946 rfield.ConOTS_local_OTS_rate[i] = 0.; 00947 } 00948 return; 00949 } 00950 00951 /* =================================================================== */ 00952 00953 /*RT_OTS_ChkSum sanity check confirms summed continua reflect contents of individuals */ 00954 void RT_OTS_ChkSum(long int ipPnt) 00955 { 00956 static long int nInsane=0; 00957 long int i; 00958 double phisig; 00959 # define LIM_INSAME_PRT 30 00960 00961 DEBUG_ENTRY( "RT_OTS_ChkSum()" ); 00962 00963 /* this entire sub is a sanity check */ 00964 /* >>chng 02 jul 23, lower bound from 0 to rfield.ipEnergyBremsThin - since now 00965 * set radiation field to black body below this energy */ 00966 for( i=rfield.ipEnergyBremsThin; i < rfield.nflux; i++ ) 00967 { 00968 phisig = rfield.otscon[i] + rfield.otslin[i] + rfield.ConInterOut[i]*rfield.lgOutOnly + 00969 rfield.outlin[i]+ 00970 rfield.outlin_noplot[i]+ 00971 rfield.ConOTS_local_OTS_rate[i]; 00972 /* >>chng 02 sep 19, add sec test on SummedDif since it can be zero whild 00973 * phisig is just above small float */ 00974 if( phisig > 0. && rfield.SummedDif[i]> 0.) 00975 { 00976 if( fabs(rfield.SummedDif[i]/phisig-1.) > 1e-3 ) 00977 { 00978 ++nInsane; 00979 /* limit amount of printout - in many failures would print entire 00980 * continuum array */ 00981 if( nInsane < LIM_INSAME_PRT ) 00982 { 00983 fprintf( ioQQQ, " PROBLEM RT_OTS_ChkSum insane SummedDif at energy %.5e error= %.2e i=%4ld\n", 00984 rfield.anu[i], rfield.SummedDif[i]/phisig - 1., i ); 00985 fprintf( ioQQQ, " SummedDif, sum are%11.4e%11.4e\n", 00986 rfield.SummedDif[i], phisig ); 00987 fprintf( ioQQQ, " otscon otslin ConInterOut outlin are%11.4e%11.4e%11.4e%11.4e\n", 00988 rfield.otscon[i], rfield.otslin[i]+rfield.outlin_noplot[i], rfield.ConInterOut[i], 00989 rfield.outlin[i]+rfield.outlin_noplot[i] ); 00990 fprintf( ioQQQ, " line continuum here are %4.4s %4.4s\n", 00991 rfield.chLineLabel[i], rfield.chContLabel[i] ); 00992 } 00993 } 00994 } 00995 00996 phisig += rfield.flux[i]; 00997 /* >>chng 02 sep 19, add sec test on SummedDif since it can be zero when 00998 * phisig is just above small float */ 00999 if( phisig > 0. && rfield.SummedDif[i]> 0. ) 01000 { 01001 if( fabs(rfield.SummedCon[i]/phisig-1.) > 1e-3 ) 01002 { 01003 ++nInsane; 01004 /* limit amount of printout - in many failures would print entire 01005 * continuum array */ 01006 if( nInsane < LIM_INSAME_PRT ) 01007 { 01008 fprintf( ioQQQ, " PROBLEM RT_OTS_ChkSum %3ld, insane SummedCon at energy %.5e error=%.2e i=%ld\n", 01009 ipPnt, rfield.anu[i], rfield.SummedCon[i]/phisig - 1., i ); 01010 fprintf( ioQQQ, " SummedCon, sum are %.4e %.4e\n", 01011 rfield.SummedCon[i], phisig ); 01012 fprintf( ioQQQ, " otscon otslin ConInterOut outlin flux are%.4e %.4e %.4e %.4e %.4e\n", 01013 rfield.otscon[i], rfield.otslin[i]+rfield.outlin_noplot[i], rfield.ConInterOut[i], 01014 rfield.outlin[i]+rfield.outlin_noplot[i], rfield.flux[i] ); 01015 fprintf( ioQQQ, " line continuum here are %s %s\n", 01016 rfield.chLineLabel[i], rfield.chContLabel[i] 01017 ); 01018 } 01019 } 01020 } 01021 } 01022 01023 if( nInsane > 0 ) 01024 { 01025 fprintf( ioQQQ, " PROBLEM RT_OTS_ChkSum too much insanity to continue.\n"); 01026 /* TotalInsanity exits after announcing a problem */ 01027 TotalInsanity(); 01028 } 01029 return; 01030 } 01031 01032 /* =================================================================== */ 01033 01034 /*RT_OTS_PrtRate print continuum and line ots rates when trace ots is on */ 01035 void RT_OTS_PrtRate( 01036 /* weakest rate to print */ 01037 double weak , 01038 /* flag, 'c' continuum, 'l' line, 'b' both */ 01039 int chFlag ) 01040 { 01041 long int i; 01042 01043 DEBUG_ENTRY( "RT_OTS_PrtRate()" ); 01044 01045 /* arg must be one of these three */ 01046 ASSERT( chFlag=='l' || chFlag=='c' || chFlag=='b' ); 01047 01048 /* 01049 * both printouts have cell number (on C array scale) 01050 * energy in ryd 01051 * the actual value of the ots rate 01052 * the ots rate relative to the continuum at that energy 01053 * rate times opacity 01054 * all are only printed if greater than weak 01055 */ 01056 01057 /*===================================================================*/ 01058 /* first print ots continua */ 01059 /*===================================================================*/ 01060 if( chFlag == 'c' || chFlag == 'b' ) 01061 { 01062 fprintf( ioQQQ, " DEBUG OTSCON array, anu, otscon, opac, OTS*opac limit:%.2e zone:%.2f IonConv?%c\n", 01063 weak,fnzone ,TorF(conv.lgConvIoniz) ); 01064 01065 for( i=0; i < rfield.nupper; i++ ) 01066 { 01067 if( rfield.otscon[i]*opac.opacity_abs[i] > weak ) 01068 { 01069 fprintf( ioQQQ, " %4ld%12.4e%12.4e%12.4e%12.4e %s \n", 01070 i, 01071 rfield.anu[i], 01072 rfield.otscon[i], 01073 opac.opacity_abs[i], 01074 rfield.otscon[i]*opac.opacity_abs[i], 01075 rfield.chContLabel[i]); 01076 01077 } 01078 } 01079 } 01080 01081 /*===================================================================*/ 01082 /* second print ots line rates */ 01083 /*===================================================================*/ 01084 if( chFlag == 'l' || chFlag == 'b' ) 01085 { 01086 fprintf( ioQQQ, " DEBUG OTSLIN array, anu, otslin, opac, OTS*opac Lab nLine limit:%.2e zone:%.2f IonConv?%c\n", 01087 weak,fnzone,TorF(conv.lgConvIoniz) ); 01088 01089 for( i=0; i < rfield.nupper; i++ ) 01090 { 01091 if( rfield.otslin[i]*opac.opacity_abs[i] > weak ) 01092 { 01093 fprintf( ioQQQ, " %4ld%12.4e%12.4e%12.4e%12.4e %s %3li\n", 01094 i, 01095 rfield.anu[i], 01096 rfield.otslin[i], 01097 opac.opacity_abs[i], 01098 rfield.otslin[i]*opac.opacity_abs[i], 01099 rfield.chLineLabel[i] , 01100 rfield.line_count[i] ); 01101 } 01102 } 01103 } 01104 return; 01105 }