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_line_driving derive radiative acceleration due to line absorption of incident continuum, 00004 * return value is line radiative acceleration */ 00005 #include "cddefines.h" 00006 #include "physconst.h" 00007 #include "iso.h" 00008 #include "dense.h" 00009 #include "taulines.h" 00010 #include "h2.h" 00011 #include "atomfeii.h" 00012 #include "rt.h" 00013 00014 /*RT_line_driving derive radiative acceleration due to line absorption of incident continuum, 00015 * return value is line radiative acceleration */ 00016 double RT_line_driving(void) 00017 { 00018 long int i, 00019 ipHi, 00020 nelem, 00021 ipLo, 00022 ipISO; 00023 00024 double AllHeavy, 00025 AllRest, 00026 OneLine, 00027 fe2drive, 00028 forlin_v, 00029 h2drive, 00030 accel_iso[NISO]; 00031 00032 /* following used for debugging */ 00033 /* double 00034 RestMax, 00035 HeavMax, 00036 hydromax; 00037 long int 00038 ipRestMax, 00039 ihmax; */ 00040 00041 DEBUG_ENTRY( "RT_line_driving()" ); 00042 00043 /* this function finds the total rate the gas absorbs energy 00044 * this result is divided by the calling routine to find the 00045 * momentum absorbed by the gas, and eventually the radiative acceleration 00046 * 00047 * the energy absorbed by the line is 00048 * Abundance * energy * A *(g_up/g_lo) * occnum * escape prob 00049 * where occnum is the photon occupation number, and the g's are 00050 * the ratios of statistical weights */ 00051 00052 /* total energy absorbed in this zone, per cubic cm 00053 * do hydrogen first */ 00054 00055 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00056 { 00057 accel_iso[ipISO] = 0; 00058 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00059 { 00060 if( (dense.IonHigh[nelem] >= nelem + 1-ipISO) ) 00061 { 00062 for( ipHi=1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ ) 00063 { 00064 /* do not put in highest level since its not real */ 00065 for( ipLo=0; ipLo < ipHi - 1; ipLo++ ) 00066 { 00067 /* do not include bogus lines */ 00068 if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont > 0 ) 00069 { 00070 OneLine = Transitions[ipISO][nelem][ipHi][ipLo].Emis->pump* 00071 Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg* 00072 Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc* 00073 dense.xIonDense[nelem][nelem+1-ipISO]; 00074 00075 accel_iso[ipISO] += OneLine; 00076 } 00077 } 00078 00079 /* do not include bogus lines */ 00080 if( SatelliteLines[ipISO][nelem][ipHi].ipCont > 0 ) 00081 { 00082 OneLine = SatelliteLines[ipISO][nelem][ipHi].Emis->pump* 00083 SatelliteLines[ipISO][nelem][ipHi].EnergyErg* 00084 SatelliteLines[ipISO][nelem][ipHi].Emis->PopOpc* 00085 dense.xIonDense[nelem][nelem+1-ipISO]; 00086 00087 accel_iso[ipISO] += OneLine; 00088 } 00089 00090 } 00091 } 00092 } 00093 accel_iso[ipISO] *= EN1RYD; 00094 } 00095 00096 /* all heavy element lines in calculation of cooling 00097 * these are the level 1 lines */ 00098 AllHeavy = 0.; 00099 for( i=1; i <= nLevel1; i++ ) 00100 { 00101 OneLine = 00102 TauLines[i].Emis->pump* 00103 TauLines[i].EnergyErg* 00104 TauLines[i].Emis->PopOpc; 00105 AllHeavy += OneLine; 00106 } 00107 00108 /* all heavy element lines treated with g-bar 00109 * these are the level 2 lines, f should be ok */ 00110 AllRest = 0.; 00111 for( i=0; i < nWindLine; i++ ) 00112 { 00113 OneLine = 00114 TauLine2[i].Emis->pump* 00115 TauLine2[i].EnergyErg* 00116 TauLine2[i].Emis->PopOpc; 00117 AllRest += OneLine; 00118 } 00119 for( i=0; i < nUTA; i++ ) 00120 { 00121 if( UTALines[i].Emis->Aul > 0. ) 00122 { 00123 OneLine = 00124 UTALines[i].Emis->pump* 00125 UTALines[i].EnergyErg* 00126 UTALines[i].Emis->PopOpc; 00127 AllRest += OneLine; 00128 } 00129 } 00130 for( i=0; i < nHFLines; i++ ) 00131 { 00132 OneLine = 00133 HFLines[i].Emis->pump* 00134 HFLines[i].EnergyErg* 00135 HFLines[i].Emis->PopOpc; 00136 AllRest += OneLine; 00137 } 00138 for( i=0; i < linesAdded2; i++) 00139 { 00140 OneLine = 00141 atmolEmis[i].pump* 00142 atmolEmis[i].tran->EnergyErg* 00143 atmolEmis[i].PopOpc; 00144 AllRest += OneLine; 00145 } 00146 for( i=0; i < nCORotate; i++ ) 00147 { 00148 OneLine = 00149 C12O16Rotate[i].Emis->pump* 00150 C12O16Rotate[i].EnergyErg* 00151 C12O16Rotate[i].Emis->PopOpc; 00152 AllRest += OneLine; 00153 OneLine = 00154 C13O16Rotate[i].Emis->pump* 00155 C13O16Rotate[i].EnergyErg* 00156 C13O16Rotate[i].Emis->PopOpc; 00157 AllRest += OneLine; 00158 } 00159 00160 /* the H2 molecule */ 00161 h2drive = H2_Accel(); 00162 00163 /* The large model FeII atom */ 00164 fe2drive = 0.; 00165 FeIIAccel(&fe2drive); 00166 00167 forlin_v = AllHeavy + accel_iso[ipH_LIKE] + accel_iso[ipHE_LIKE] + 00168 fe2drive + h2drive + AllRest; 00169 00170 /*fprintf(ioQQQ," wind te %e %e %e %e %e\n", 00171 AllHeavy , HydroAccel , fe2drive , he1l , AllRest );*/ 00172 return( forlin_v ); 00173 }