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 /*ParseMagnet parse magnetic field command */ 00004 /*Magnetic_init initialize magnetic field parameters */ 00005 /*Magnetic_reinit - reinitialized magnetic field at start of new iteration */ 00006 /*Magnetic_evaluate evaluate some parameters to do with magnetic field */ 00007 #include "cddefines.h" 00008 #include "physconst.h" 00009 #include "dense.h" 00010 #include "doppvel.h" 00011 #include "optimize.h" 00012 #include "input.h" 00013 #include "wind.h" 00014 #include "magnetic.h" 00015 00016 /* the initial magnetic field */ 00017 static double Btangl_init; 00018 00019 /* this is logical var, set in zero, which says whether the magnetic field has been 00020 * initialized */ 00021 static bool lgBinitialized; 00022 00023 /* the current magnetic field */ 00024 static double Btangl_here; 00025 00026 /* the initial parallel and tangential fields for ordered case */ 00027 static double Bpar_init, Btan_init; 00028 00029 /* the current parallel and tangential fields for ordered case */ 00030 static double Bpar_here, Btan_here; 00031 00032 /* this is the gamma power law index, default is 4. / 3. */ 00033 static double gamma_mag; 00034 00035 /*Magnetic_evaluate evaluate some parameters to do with magnetic field */ 00036 void Magnetic_evaluate(void) 00037 { 00038 00039 DEBUG_ENTRY( "Magnetic_evaluate()" ); 00040 00041 /* this flag set true if magnetic field is specified */ 00042 if( magnetic.lgB ) 00043 { 00044 static double density_initial, 00045 /* the square of the Alfven velocity at illuminated face */ 00046 v_A; 00047 00048 /* does magnetic field need to be initialized for this iteration? 00049 * flag is reset false at init of code, and at start of every iteration */ 00050 if( !lgBinitialized ) 00051 { 00052 lgBinitialized = true; 00053 00054 /* set initial tangled field */ 00055 Btangl_here = Btangl_init; 00056 00057 /* set initial ordered field */ 00058 /* mag field angle_wrt_los set when ordered field specified */ 00059 Bpar_here = Bpar_init; 00060 Btan_here = Btan_init; 00061 00062 /* XMassDensity was set above, safe to use this on first call */ 00063 density_initial = dense.xMassDensity; 00064 00065 /* this is used to update tangential field */ 00066 v_A = POW2(Bpar_init) / (PI4 * density_initial ); 00067 } 00068 00069 /* now update parameters in tangled field case */ 00070 /* magnetic pressure is a function of the local density, use gamma law */ 00071 Btangl_here = Btangl_init * pow(dense.xMassDensity/density_initial, gamma_mag/2.); 00072 00073 /* ordered components of field - parallel field is always constant - find tangential component - 00074 * but only in wind case */ 00075 if( wind.windv0 != 0. ) 00076 { 00077 /* N B - must preserve sign in this equation - will blow if product of wind speeds is equal to v_A) */ 00078 /* wind.windv*wind.windv0 == v_A should not occur since mag pressure goes to inf */ 00079 Btan_here = Btan_init * (POW2(wind.windv0) - v_A)/ (wind.windv*wind.windv0-v_A); 00080 } 00081 00082 /* magnetic pressure - sum of two fields - can have negative pressure (tension) 00083 * is ordered field dominates */ 00084 magnetic.pressure = POW2(Btangl_here)/PI8 + (POW2(Btan_here) - POW2(Bpar_here)) / PI8; 00085 00086 /* energy density - this is positive */ 00087 magnetic.energydensity = POW2(Btangl_here)/PI8 + (POW2(Btan_here) + POW2(Bpar_here)) / PI8; 00088 00089 /* option for turbulence in equipartition with B field */ 00090 if( DoppVel.lgTurbEquiMag ) 00091 { 00092 /* >>chng 05 jan 26, as per Robin Williams email, 00093 * evaluate energydensity above, which is +ve, and use that for 00094 * velocity here - had used pressure but could not evaluate when negative */ 00095 /* turbulent velocity is mag pres over density */ 00096 /* >>chng 06 apr 19, use DoppVel.Heiles_Troland_F */ 00097 /*DoppVel.TurbVel = (realnum)sqrt(2.*magnetic.energydensity/dense.xMassDensity);*/ 00098 DoppVel.TurbVel = (realnum)sqrt(6.*magnetic.energydensity/dense.xMassDensity/ 00099 DoppVel.Heiles_Troland_F); 00100 /* >>chng 06 apr 19, do not double mag pressure, really count turbulence as pressure */ 00101 /* double magnetic pressure to account for ram pressure due to turbulence, 00102 * which is not counted elsewhere 00103 magnetic.pressure *= 2.;*/ 00104 } 00105 00106 /* input parser made sure gamma != 1, default magnetic gamma is 4/3 */ 00107 magnetic.EnthalpyDensity = gamma_mag/(gamma_mag-1.) * 00108 POW2(Btangl_here)/PI8 + (POW2(Btan_here) + POW2(Bpar_here)) / PI4; 00109 } 00110 else 00111 { 00112 magnetic.pressure = 0.; 00113 magnetic.energydensity = 0.; 00114 magnetic.EnthalpyDensity = 0.; 00115 } 00116 return; 00117 } 00118 00119 /*Magnetic_reinit - reinitialized magnetic field at start of new iteration */ 00120 void Magnetic_reinit(void) 00121 { 00122 DEBUG_ENTRY( "Magnetic_reinit()" ); 00123 00124 /* this says whether B has been initialized in this run */ 00125 lgBinitialized = false; 00126 return; 00127 } 00128 00129 /* Magnetic_init initialize magnetic field parameters */ 00130 void Magnetic_init(void) 00131 { 00132 00133 DEBUG_ENTRY( "Magnetic_init()" ); 00134 00135 gamma_mag = 4./3.; 00136 magnetic.lgB = false; 00137 /* this says whether B has been initialized in this run */ 00138 lgBinitialized = false; 00139 /* the initial tangled and ordered fields */ 00140 Btangl_init = 0.; 00141 Btangl_here = DBL_MAX; 00142 magnetic.pressure = DBL_MAX; 00143 magnetic.energydensity = DBL_MAX; 00144 Bpar_init = 0.; 00145 Btan_init = 0.; 00146 Bpar_here = DBL_MAX; 00147 Btan_here = DBL_MAX; 00148 magnetic.EnthalpyDensity = DBL_MAX; 00149 return; 00150 } 00151 00152 /*ParseMagnet parse magnetic field command */ 00153 void ParseMagnet(char *chCard ) 00154 { 00155 bool lgEOL; 00156 long int i; 00157 bool lgTangled; 00158 double angle_wrt_los=-1. , Border_init=-1.; 00159 00160 DEBUG_ENTRY( "ParseMagnet()" ); 00161 00162 /* flag saying B turned on */ 00163 magnetic.lgB = true; 00164 00165 /* check whether ordered is specified - if not this is tangled */ 00166 if( nMatch("ORDE" , chCard ) ) 00167 { 00168 /* ordered field case */ 00169 lgTangled = false; 00170 00171 /* field is ordered, not isotropic - need field direction - spcified 00172 * by angle away from beam of light - 0 => parallel to beam */ 00173 00174 i = 5; 00175 /* this is the log of the ordered field strength */ 00176 Border_init = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00177 Border_init = pow(10.,Border_init ); 00178 00179 /* this is the angle (default in degrees) with respect to the line of sight */ 00180 angle_wrt_los = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00181 if( lgEOL ) 00182 NoNumb(chCard); 00183 00184 /* if radians is on the line then the number already is in radians - 00185 * else convert to radians */ 00186 if( !nMatch("RADI" , chCard ) ) 00187 angle_wrt_los *= PI2 / 360.; 00188 00189 /* now get initial components that we will work with */ 00190 Bpar_init = Border_init*cos(angle_wrt_los); 00191 Btan_init = Border_init*sin(angle_wrt_los); 00192 00193 } 00194 else 00195 { 00196 /* tangled field case */ 00197 lgTangled = true; 00198 i = 5; 00199 /* this is the log of the tangled field strength */ 00200 Btangl_init = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00201 if( lgEOL ) 00202 NoNumb(chCard); 00203 Btangl_init = pow(10.,Btangl_init ); 00204 00205 /* optional gamma for dependence on pressure */ 00206 gamma_mag = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00207 if( lgEOL ) 00208 gamma_mag = 4./3.; 00209 00210 if( gamma_mag !=0. && gamma_mag <=1. ) 00211 { 00212 /* impossible value for gamma */ 00213 fprintf( ioQQQ, 00214 " This value of gamma (%.3e) is impossible. Must have gamma = 0 or > 1.\n Sorry.\n", 00215 gamma_mag ); 00216 cdEXIT(EXIT_FAILURE); 00217 } 00218 } 00219 00220 /* vary option */ 00221 if( optimize.lgVarOn ) 00222 { 00223 /* number of parameters */ 00224 optimize.nvarxt[optimize.nparm] = 2; 00225 if( lgTangled ) 00226 { 00227 /* tangled field case */ 00228 strcpy( optimize.chVarFmt[optimize.nparm], "MAGNETIC FIELD TANGLED =%f GAMMA= %f" ); 00229 optimize.vparm[0][optimize.nparm] = (realnum)log10( Btangl_init ); 00230 /* this is not varied */ 00231 optimize.vparm[1][optimize.nparm] = (realnum)gamma_mag; 00232 } 00233 else 00234 { 00235 /* ordered field case */ 00236 strcpy( optimize.chVarFmt[optimize.nparm], "MAGNETIC FIELD ORDERED =%f ANGLE RADIANS = %f" ); 00237 optimize.vparm[0][optimize.nparm] = (realnum)log10( Border_init ); 00238 /* this is not varied */ 00239 optimize.vparm[1][optimize.nparm] = (realnum)angle_wrt_los; 00240 } 00241 00242 /* pointer to where to write */ 00243 optimize.nvfpnt[optimize.nparm] = input.nRead; 00244 optimize.vincr[optimize.nparm] = 0.5; 00245 ++optimize.nparm; 00246 } 00247 return; 00248 }