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