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