00001
00002
00003
00004
00005
00006
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 #include "parser.h"
00016
00017 t_magnetic magnetic;
00018
00019
00020 static double Btangl_init;
00021
00022
00023
00024 static bool lgBinitialized;
00025
00026
00027 static double Btangl_here;
00028
00029
00030 static double Bpar_init, Btan_init;
00031
00032
00033 static double Bpar_here, Btan_here;
00034
00035
00036 static double gamma_mag;
00037
00038
00039 void Magnetic_evaluate(void)
00040 {
00041
00042 DEBUG_ENTRY( "Magnetic_evaluate()" );
00043
00044
00045 if( magnetic.lgB )
00046 {
00047 static double density_initial,
00048
00049 v_A;
00050
00051
00052
00053 if( !lgBinitialized )
00054 {
00055 lgBinitialized = true;
00056
00057
00058 Btangl_here = Btangl_init;
00059
00060
00061
00062 Bpar_here = Bpar_init;
00063 Btan_here = Btan_init;
00064
00065
00066 density_initial = dense.xMassDensity;
00067
00068
00069 v_A = POW2(Bpar_init) / (PI4 * density_initial );
00070 }
00071
00072
00073
00074 Btangl_here = Btangl_init * pow(dense.xMassDensity/density_initial, gamma_mag/2.);
00075
00076
00077
00078 if( !wind.lgStatic() )
00079 {
00080
00081
00082 Btan_here = Btan_init * (POW2(wind.windv0) - v_A)/ (wind.windv*wind.windv0-v_A);
00083 }
00084
00085
00086
00087 magnetic.pressure = POW2(Btangl_here)/PI8 + (POW2(Btan_here) - POW2(Bpar_here)) / PI8;
00088
00089
00090 magnetic.energydensity = POW2(Btangl_here)/PI8 + (POW2(Btan_here) + POW2(Bpar_here)) / PI8;
00091
00092
00093 if( DoppVel.lgTurbEquiMag )
00094 {
00095
00096
00097
00098
00099
00100
00101 DoppVel.TurbVel = (realnum)sqrt(6.*magnetic.energydensity/dense.xMassDensity/
00102 DoppVel.Heiles_Troland_F);
00103
00104
00105
00106
00107 }
00108
00109
00110 magnetic.EnthalpyDensity = gamma_mag/(gamma_mag-1.) *
00111 POW2(Btangl_here)/PI8 + (POW2(Btan_here) + POW2(Bpar_here)) / PI4;
00112 }
00113 else
00114 {
00115 magnetic.pressure = 0.;
00116 magnetic.energydensity = 0.;
00117 magnetic.EnthalpyDensity = 0.;
00118 }
00119 return;
00120 }
00121
00122
00123 void Magnetic_reinit(void)
00124 {
00125 DEBUG_ENTRY( "Magnetic_reinit()" );
00126
00127
00128 lgBinitialized = false;
00129 return;
00130 }
00131
00132
00133 void Magnetic_init(void)
00134 {
00135
00136 DEBUG_ENTRY( "Magnetic_init()" );
00137
00138 gamma_mag = 4./3.;
00139 magnetic.lgB = false;
00140
00141 lgBinitialized = false;
00142
00143 Btangl_init = 0.;
00144 Btangl_here = DBL_MAX;
00145 magnetic.pressure = DBL_MAX;
00146 magnetic.energydensity = DBL_MAX;
00147 Bpar_init = 0.;
00148 Btan_init = 0.;
00149 Bpar_here = DBL_MAX;
00150 Btan_here = DBL_MAX;
00151 magnetic.EnthalpyDensity = DBL_MAX;
00152 return;
00153 }
00154
00155
00156 void ParseMagnet(Parser &p )
00157 {
00158 bool lgTangled;
00159 double angle_wrt_los=-1. , Border_init=-1.;
00160
00161 DEBUG_ENTRY( "ParseMagnet()" );
00162
00163
00164 magnetic.lgB = true;
00165
00166
00167 if( p.nMatch("ORDE") )
00168 {
00169
00170 lgTangled = false;
00171
00172
00173
00174
00175
00176 Border_init = p.getNumberCheckAlwaysLog("ordered field");
00177
00178
00179 angle_wrt_los = p.getNumberCheck("LOS angle");
00180
00181
00182
00183 if( !p.nMatch("RADI") )
00184 angle_wrt_los *= PI2 / 360.;
00185
00186
00187 Bpar_init = Border_init*cos(angle_wrt_los);
00188 Btan_init = Border_init*sin(angle_wrt_los);
00189
00190 }
00191 else
00192 {
00193
00194 lgTangled = true;
00195
00196 Btangl_init = p.getNumberCheckAlwaysLog("tangled field");
00197
00198
00199 gamma_mag = p.getNumberDefault("field gamma law", 4./3.);
00200
00201 if( gamma_mag != 0. && gamma_mag <= 1. )
00202 {
00203
00204 fprintf( ioQQQ,
00205 " This value of gamma (%.3e) is impossible. Must have gamma = 0 or > 1.\n Sorry.\n",
00206 gamma_mag );
00207 cdEXIT(EXIT_FAILURE);
00208 }
00209 }
00210
00211
00212 if( optimize.lgVarOn )
00213 {
00214
00215 optimize.nvarxt[optimize.nparm] = 2;
00216 if( lgTangled )
00217 {
00218
00219
00220 strcpy( optimize.chVarFmt[optimize.nparm], "MAGNETIC FIELD TANGLED= %f LOG GAMMA= %f" );
00221 optimize.vparm[0][optimize.nparm] = (realnum)log10( Btangl_init );
00222
00223 optimize.vparm[1][optimize.nparm] = (realnum)gamma_mag;
00224 }
00225 else
00226 {
00227
00228
00229 strcpy( optimize.chVarFmt[optimize.nparm], "MAGNETIC FIELD ORDERED= %f LOG ANGLE RADIANS= %f" );
00230 optimize.vparm[0][optimize.nparm] = (realnum)log10( Border_init );
00231
00232 optimize.vparm[1][optimize.nparm] = (realnum)angle_wrt_los;
00233 }
00234
00235
00236 optimize.nvfpnt[optimize.nparm] = input.nRead;
00237 optimize.vincr[optimize.nparm] = 0.5;
00238 ++optimize.nparm;
00239 }
00240 return;
00241 }