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 /*PrtLinePres print line radiation pressures for current conditions */ 00004 #include "cddefines.h" 00005 #include "pressure.h" 00006 #include "taulines.h" 00007 #include "iso.h" 00008 #include "dense.h" 00009 #include "lines_service.h" 00010 #include "h2.h" 00011 #include "prt.h" 00012 /* faintest pressure we will bother with */ 00013 #define THRESH 0.05 00014 00015 void PrtLinePres(void) 00016 { 00017 long int i, 00018 ip, 00019 ipLo, 00020 ipHi, 00021 nelem; 00022 int ier; 00023 double RadPres1; 00024 00025 /* this is limit to number of lines we can store at one time */ 00026 # define NLINE 100 00027 /* labels, wavelengths, and fraction of total pressure, for important lines */ 00028 char chLab[NLINE][5]; 00029 realnum wl[NLINE] , frac[NLINE]; 00030 long int iperm[NLINE]; 00031 00032 /* will be used to check on size of opacity, was capped at this value */ 00033 realnum smallfloat=SMALLFLOAT*10.f; 00034 00035 DEBUG_ENTRY( "PrtLinePres()" ); 00036 00037 /* this routine only called if printout of contributors to line 00038 * radiation pressure is desired */ 00039 00040 /* compute radiation pressure due to iso lines */ 00041 ip = 0; 00042 00043 for(i=0; i<NLINE; ++i) 00044 { 00045 frac[i] = FLT_MAX; 00046 wl[i] = FLT_MAX; 00047 } 00048 00049 if( pressure.pres_radiation_lines_curr > 1e-30 ) 00050 { 00054 for( long ipISO = 0; ipISO<NISO; ipISO++ ) 00055 { 00056 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00057 { 00058 /* does this ion stage exist? */ 00059 if( dense.IonHigh[nelem] >= nelem + 1 - ipISO ) 00060 { 00061 /* do not include highest levels since maser can occur due to topoff, 00062 * and pops are set to small number in this case */ 00063 for( ipHi=1; ipHi <iso.numLevels_local[ipISO][nelem] - iso.nCollapsed_local[ipISO][nelem]; ipHi++ ) 00064 { 00065 for( ipLo=0; ipLo < ipHi; ipLo++ ) 00066 { 00067 if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 ) 00068 continue; 00069 00070 ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul > iso.SmallA ); 00071 00072 /*NB - this code must be kept parallel with code in pressure_total */ 00073 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc > smallfloat && 00074 /* >>chng 01 nov 01, add test that have not overrun optical depth scale */ 00075 ( (Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauTot - Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn) > smallfloat ) ) 00076 { 00077 RadPres1 = PressureRadiationLine( &Transitions[ipISO][nelem][ipHi][ipLo], dense.xIonDense[nelem][nelem+1-ipISO] ); 00078 00079 if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH ) 00080 { 00081 wl[ip] = Transitions[ipISO][nelem][ipHi][ipLo].WLAng; 00082 00083 /* put null terminated line label into chLab using line array*/ 00084 chIonLbl(chLab[ip], &Transitions[ipISO][nelem][ipHi][ipLo]); 00085 frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr); 00086 00087 ip = MIN2((long)NLINE-1,ip+1); 00088 } 00089 } 00090 } 00091 } 00092 } 00093 } 00094 } 00095 00096 /* line radiation pressure from large set of level 1 lines */ 00097 for( i=1; i <= nLevel1; i++ ) 00098 { 00099 if( TauLines[i].Hi->Pop > 1e-30 ) 00100 { 00101 RadPres1 = PressureRadiationLine( &TauLines[i], 1. ); 00102 } 00103 else 00104 { 00105 RadPres1 = 0.; 00106 } 00107 00108 if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH ) 00109 { 00110 wl[ip] = TauLines[i].WLAng; 00111 00112 /* put null terminated line label into chLab using line array*/ 00113 chIonLbl(chLab[ip], &TauLines[i]); 00114 frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr); 00115 00116 ip = MIN2((long)NLINE-1,ip+1); 00117 } 00118 } 00119 00120 for( i=0; i < nWindLine; i++ ) 00121 { 00122 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO ) 00123 { 00124 if( TauLine2[i].Hi->Pop > 1e-30 ) 00125 { 00126 RadPres1 = PressureRadiationLine( &TauLine2[i], 1. ); 00127 } 00128 else 00129 { 00130 RadPres1 = 0.; 00131 } 00132 00133 if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH ) 00134 { 00135 wl[ip] = TauLine2[i].WLAng; 00136 00137 /* put null terminated line label into chLab using line array*/ 00138 chIonLbl(chLab[ip], &TauLine2[i]); 00139 frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr); 00140 00141 ip = MIN2((long)NLINE-1,ip+1); 00142 } 00143 } 00144 } 00145 00146 for( i=0; i < nHFLines; i++ ) 00147 { 00148 if( HFLines[i].Hi->Pop > 1e-30 ) 00149 { 00150 RadPres1 = PressureRadiationLine( &HFLines[i], 1. ); 00151 } 00152 else 00153 { 00154 RadPres1 = 0.; 00155 } 00156 00157 if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH ) 00158 { 00159 wl[ip] = HFLines[i].WLAng; 00160 00161 /* put null terminated line label into chLab using line array*/ 00162 chIonLbl(chLab[ip], &HFLines[i]); 00163 frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr); 00164 00165 ip = MIN2((long)NLINE-1,ip+1); 00166 } 00167 } 00168 00169 /*Chiani & Leiden Atoms & Molecules*/ 00170 for( i=0; i <linesAdded2; i++ ) 00171 { 00172 if( atmolEmis[i].tran->Hi->Pop > 1e-30 ) 00173 { 00174 RadPres1 = PressureRadiationLine( atmolEmis[i].tran, 1. ); 00175 } 00176 else 00177 { 00178 RadPres1 = 0.; 00179 } 00180 00181 if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH ) 00182 { 00183 wl[ip] = atmolEmis[i].tran->WLAng; 00184 00185 /* put null terminated line label into chLab using line array*/ 00186 chIonLbl(chLab[ip], atmolEmis[i].tran); 00187 frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr); 00188 00189 ip = MIN2((long)NLINE-1,ip+1); 00190 } 00191 } 00192 00193 /* radiation pressure due to H2 */ 00194 RadPres1 = H2_RadPress(); 00195 if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH ) 00196 { 00197 wl[ip] = 0; 00198 00199 /* put null terminated 4 char line label into chLab using line array*/ 00200 strcpy(chLab[ip], "H2 "); 00201 frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr); 00202 00203 ip = MIN2((long)NLINE-1,ip+1); 00204 } 00205 00206 for( i=0; i < nCORotate; i++ ) 00207 { 00208 if( C12O16Rotate[i].Hi->Pop > 1e-30 ) 00209 { 00210 RadPres1 = PressureRadiationLine( &C12O16Rotate[i], 1. ); 00211 } 00212 else 00213 { 00214 RadPres1 = 0.; 00215 } 00216 00217 if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH ) 00218 { 00219 wl[ip] = C12O16Rotate[i].WLAng; 00220 00221 /* put null terminated line label into chLab using line array*/ 00222 chIonLbl(chLab[ip], &C12O16Rotate[i]); 00223 frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr); 00224 00225 ip = MIN2((long)NLINE-1,ip+1); 00226 } 00227 } 00228 00229 for( i=0; i < nCORotate; i++ ) 00230 { 00231 if( C13O16Rotate[i].Hi->Pop > 1e-30 ) 00232 { 00233 RadPres1 = PressureRadiationLine( &C13O16Rotate[i], 1. ); 00234 } 00235 else 00236 { 00237 RadPres1 = 0.; 00238 } 00239 00240 if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH ) 00241 { 00242 wl[ip] = C13O16Rotate[i].WLAng; 00243 00244 /* put null terminated line label into chLab using line array*/ 00245 chIonLbl(chLab[ip], &C13O16Rotate[i]); 00246 frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr); 00247 00248 ip = MIN2((long)NLINE-1,ip+1); 00249 } 00250 } 00251 00252 /* return if no significant contributors to radiation pressure found */ 00253 if( ip<= 0 ) 00254 return; 00255 00256 /* sort from strong to weak, then only print strongest */ 00257 spsort( 00258 /* input array to be sorted */ 00259 frac, 00260 /* number of values in x */ 00261 ip, 00262 /* permutation output array */ 00263 iperm, 00264 /* flag saying what to do - 1 sorts into increasing order, not changing 00265 * the original routine */ 00266 -1, 00267 /* error condition, should be 0 */ 00268 &ier); 00269 00270 /* now print up to ten strongest contributors to radiation pressure */ 00271 fprintf( ioQQQ, " P(Lines):" ); 00272 for( i=0; i < MIN2(10,ip); i++ ) 00273 { 00274 int ipline = iperm[i]; 00275 fprintf( ioQQQ, "(%4.4s ", chLab[ipline]); 00276 prt_wl(ioQQQ, wl[ipline] ); 00277 fprintf(ioQQQ," %.2f) ",frac[ipline]); 00278 } 00279 00280 /* finally end the line */ 00281 fprintf( ioQQQ, "\n" ); 00282 } 00283 return; 00284 }