00001
00009 #include "DiffusionCoefficient.h"
00010
00011 #include <math.h>
00012 #include <vector>
00013 #include "rroots.h"
00014 #include "../VariousFunctions/bisection.h"
00015
00016 #include "../Logging/Output.h"
00017 #include "../Exceptions/error.h"
00018
00019 #include <ctime>
00020
00021 using namespace std;
00022
00023 #define double_zero 1.e-21
00024 #define min_Dxx 1.e-21
00025
00026 #if defined(_WIN32) || defined(_WIN64)
00027 #define strncasecmp _strnicmp
00028 #define strcasecmp _stricmp
00029 #endif
00030
00031
00037 DiffusionCoefficient::DiffusionCoefficient(Grid &grid, DiffusionCoefficientParamStructure DxxParamStructure) {
00038
00039 this->change_ind = clock();
00040
00041 this->AllocateMemory(grid.L.size, grid.pc.size, grid.alpha.size);
00042
00043 this->Get(grid, DxxParamStructure);
00044 }
00045
00050 DiffusionCoefficient::DiffusionCoefficient(Grid &grid) : Matrix3D<double>(grid.L.size, grid.pc.size, grid.alpha.size) {
00051
00052 this->change_ind = clock();
00053
00054 initialized = false;
00055 }
00056
00062 void DiffusionCoefficient::AllocateMemory(int L_size, int pc_size, int alpha_size) {
00063
00064
00065 Matrix3D<double>::AllocateMemory(L_size, pc_size, alpha_size);
00066 }
00067
00068
00076 void DiffusionCoefficient::Get(Grid &grid, DiffusionCoefficientParamStructure DxxParamStructure) {
00077 if (!Matrix3D<double>::initialized) {
00078 throw error_msg("DXX_UNINITIALIZED_MATRIX", "Matrix has not been initialized.");
00079 }
00080 this->name = DxxParamStructure.DxxName + " " + DxxParamStructure.waveName + "";
00081 Dxx_parameters = DxxParamStructure;
00082
00083 this->type = DxxParamStructure.DxxType;
00084 this->useScale = DxxParamStructure.useScale;
00085
00086
00087
00088 if (DxxParamStructure.loadOrCalculate == "LOC_LOAD" || DxxParamStructure.loadOrCalculate == "LOC_LOADORCALCULATE") {
00089 try {
00090
00091 LoadDiffusionCoefficient(grid.L, grid.epc, grid.alpha, DxxParamStructure.filename.c_str(), DxxParamStructure.filetype);
00092 } catch (error_msg err) {
00093 if (DxxParamStructure.loadOrCalculate == "LOC_LOADORCALCULATE") {
00094 Output::echo(2, "%s", err.what().c_str());
00095
00096 Calculate(grid.L, grid.epc, grid.alpha, DxxParamStructure);
00097
00098 (*this) = (*this) * (60.0 * 60 * 24 * DxxParamStructure.MLT_averaging);
00099
00100 this->writeToFile(DxxParamStructure.filename.c_str(), grid.L, grid.epc, grid.alpha);
00101 } else {
00102 err.add("UNKNOWN_ERROR", "Unknown error while loading diffusion coefficient file");
00103 throw err;
00104 }
00105 }
00106 } else if (DxxParamStructure.loadOrCalculate == "LOC_CALCULATE" || DxxParamStructure.loadOrCalculate == "LOC_LOADORCALCULATE") {
00107 Calculate(grid.L, grid.epc, grid.alpha, DxxParamStructure);
00108
00109 (*this) = (*this) * (60.0 * 60 * 24 * DxxParamStructure.MLT_averaging);
00110
00111 this->writeToFile(DxxParamStructure.filename.c_str(), grid.L, grid.epc, grid.alpha);
00112 } else {
00113 throw error_msg("DXX_GET", "Undefined load or calculate choise.");
00114 }
00115
00116
00117 if (DxxParamStructure.multiplicator != 1) {
00118 (*this) = (*this) * DxxParamStructure.multiplicator;
00119 Output::echo(3, "%s is multiplied by %f.\n", this->name.c_str(), DxxParamStructure.multiplicator);
00120 }
00121 is_active = false;
00122
00123
00124 if (this->type == "DCT_Dpp" || this->type == "DCT_DppLpp") (*this) = (*this) * grid.pc * grid.pc;
00125 if (this->type == "DCT_Dpa" || this->type == "DCT_DpaLpp") (*this) = (*this) * grid.pc;
00126
00127 initialized = true;
00128 }
00129
00136 void DiffusionCoefficient::Calculate(GridElement &L, GridElement &epc, GridElement &Alpha, DiffusionCoefficientParamStructure DxxParamStructure) {
00137 int il, im, ia;
00138 bool positivonly=true;
00139 double (*int_Dxx_loc)(double lambda, DiffusionCoefficientParamStructure DxxParamStructure);
00140
00141
00142 if (DxxParamStructure.DxxType == "DCT_Daa" || DxxParamStructure.DxxType == "DCT_DaaLpp") {
00143 int_Dxx_loc = int_Daa_loc;
00144 } else if (DxxParamStructure.DxxType == "DCT_Dpp" || DxxParamStructure.DxxType == "DCT_DppLpp") {
00145 int_Dxx_loc = int_Dpp_loc;
00146 } else if (DxxParamStructure.DxxType == "DCT_Dpa" || DxxParamStructure.DxxType == "DCT_DpaLpp") {
00147 int_Dxx_loc = int_Dpa_loc;
00148 positivonly = false;
00149 } else {
00150 throw error_msg("DXX_TYPE", "Unknown diffusion coefficient");
00151 }
00152
00153 Output::echo(2, "Calculating %s for %s... ", DxxParamStructure.DxxName.c_str(), DxxParamStructure.waveName.c_str());
00154
00155 for (il = 0; il < L.size ; il++) {
00156 for (im = 0; im < epc.size; im++) {
00157 for (ia = 0; ia < Alpha.size; ia++) {
00158
00159
00160 (*this)[il][im][ia] = Dxx_ba(L[il][im][ia], epc[il][im][ia]/VC::mc2, Alpha[il][im][ia], int_Dxx_loc, DxxParamStructure) / 2;
00161
00162 if (positivonly && ((*this)[il][im][ia] < min_Dxx)) {
00163 (*this)[il][im][ia] = min_Dxx;
00164 } else if (!((*this)[il][im][ia] >= 0) && !((*this)[il][im][ia] <= 0)) {
00165 if (positivonly) {
00166 (*this)[il][im][ia] = min_Dxx;
00167 } else {
00168 (*this)[il][im][ia] = 0;
00169 }
00170 }
00171 }
00172 }
00173 }
00174
00175 Output::echo(2, "done!\n");
00176 }
00177
00178
00182 bool DiffusionCoefficient::LoadDiffusionCoefficient(GridElement &L, GridElement &pc, GridElement &alpha, string D_filename, string filetype) {
00183 if (!Matrix3D<double>::initialized) {
00184 throw error_msg("DXX_UNINITIALIZED_MATRIX", "Matrix has not been initialized.");
00185 }
00186 Output::echo(2, "Loading %s... ", D_filename.c_str());
00187
00188
00189 if (filetype == "IFT_GRID") {
00190 LoadDiffusionCoefficientFromFileWithGrid(L, pc, alpha, D_filename);
00191 } else if (filetype == "IFT_GRID_LPA") {
00192 LoadDiffusionCoefficientFromFileWithGrid(L, pc, alpha, D_filename, "IFT_GRID_LPA");
00193 } else if (filetype == "IFT_GRID_APL") {
00194 LoadDiffusionCoefficientFromFileWithGrid(L, pc, alpha, D_filename, "IFT_GRID_APL");
00195 } else if (filetype == "IFT_PLANE") {
00196
00197 LoadDiffusionCoefficientFromPlaneFile(L, pc, alpha, D_filename);
00198 } else {
00199 throw error_msg("DXX_LOAD_TYPE_ERR", "Unknown coefficients file format %s", filetype.c_str());
00200 }
00201 Output::echo(2, "done!\n");
00202 return true;
00203 }
00204
00207 bool DiffusionCoefficient::LoadDiffusionCoefficientFromPlaneFile(GridElement &L, GridElement &pc, GridElement &alpha, string D_filename) {
00208
00209 int il, im, ia;
00210
00211
00212 ifstream input_D(D_filename.c_str());
00213
00214 if (input_D != NULL && !input_D.eof()) {
00215 for (il = 0; il < L.size; il++)
00216 for (im = 0; im < pc.size; im++)
00217 for (ia = 0; ia < alpha.size; ia++) {
00218
00219 input_D >> (*this)[il][im][ia];
00220 if (input_D == NULL) {
00221 throw error_msg("LOAD_ERROR", "Wrong coefficient file %s (to few values)\n", D_filename.c_str());
00222 }
00223 }
00224 } else {
00225 throw error_msg("DXX_LOAD_GRID_ERR", "error reading input file %s\n", D_filename.c_str());
00226
00227 return false;
00228 }
00229 input_D.close();
00230 return true;
00231
00232 }
00233
00234
00237 bool DiffusionCoefficient::LoadDiffusionCoefficientFromFileWithGrid(GridElement &L, GridElement &epc, GridElement &alpha, string D_filename, string gridOrder) {
00238 int il, im, ia;
00239 double Load_L, Load_epc, Load_alpha;
00240 double err = 1.e-3;
00241 string tmp;
00242
00243
00244 ifstream input_D(D_filename.c_str());
00245
00246 if (input_D == NULL) {
00247 throw error_msg("DXX_FILE_OPEN_ERR", "Diffusion coefficient file %s not found.", D_filename.c_str());
00248 }
00249
00250
00251 input_D >> tmp;
00252 while (strcasecmp(tmp.c_str(), "ZONE") != 0 && !input_D.eof() ) {
00253 input_D >> tmp;
00254 Output::echo(2, "%s ", tmp.c_str());
00255 if (input_D.eof()) {
00256 throw error_msg("DXX_LOAD_GRID_ERR", "Keyword 'ZONE' not found in the diffusion coefficient file: %s\n", D_filename.c_str());
00257 }
00258 }
00259
00260 input_D.ignore(9999, '\n');
00261
00262 if (!input_D.eof()) {
00263
00264 if (gridOrder == "IFT_GRID_LPA") {
00265 for (il = 0; il < L.size; il++)
00266 for (im = 0; im < epc.size; im++)
00267 for (ia = 0; ia < alpha.size; ia++) {
00268 input_D >> Load_L >> Load_epc >> Load_alpha;
00269
00270
00271 input_D >> (*this)[il][im][ia];
00272
00273 if (fabs(Load_L - L[il][im][ia]) > err || fabs(Load_epc - epc[il][im][ia]) > err || fabs(Load_alpha - alpha[il][im][ia]) > err) {
00274 throw error_msg("DXX_LOAD_GRID_ERR", "Loading %s: grid mismatch.\nLoaded: L=%e, epc=%e, alpha=%e\nGrid: L=%e, epc=%e, alpha=%e\n", D_filename.c_str(), Load_L, Load_epc, Load_alpha, L[il][im][ia], epc[il][im][ia], alpha[il][im][ia]);
00275 }
00276 }
00277 } else if (gridOrder == "IFT_GRID_APL") {
00278 for (ia = 0; ia < alpha.size; ia++)
00279 for (im = 0; im < epc.size; im++)
00280 for (il = 0; il < L.size; il++) {
00281 input_D >> Load_alpha >> Load_epc >> Load_L ;
00282
00283
00284 input_D >> (*this)[il][im][ia];
00285
00286 if (fabs(Load_L - L[il][im][ia]) > err || fabs(Load_epc - epc[il][im][ia]) > err || fabs(Load_alpha - alpha[il][im][ia]) > err) {
00287 throw error_msg("DXX_LOAD_GRID_ERR", "Loading %s: grid mismatch.\nLoaded: L=%e, epc=%e, alpha=%e\nGrid: L=%e, epc=%e, alpha=%e\n", D_filename.c_str(), Load_L, Load_epc, Load_alpha, L[il][im][ia], epc[il][im][ia], alpha[il][im][ia]);
00288 }
00289 }
00290 } else {
00291 throw error_msg("UNKNOWN_GRID_TYPE", "Unknown grid order: %s", gridOrder.c_str());
00292 }
00293 }
00294 input_D.close();
00295 return true;
00296
00297 }
00298
00301 void DiffusionCoefficient::MakeDLL(GridElement &L, GridElement &pc, GridElement &alpha, double Kp, string DLLType) {
00302 if (DLLType == "DLLT_BE") {
00303 MakeDLL_BE(L, pc, alpha, Kp);
00304 } else if (DLLType == "DLLT_B") {
00305 MakeDLL_B(L, pc, alpha, Kp);
00306 } else if (DLLType == "DLLT_FAKE") {
00307 MakeDLL_FAKE(L, pc, alpha, Kp);
00308 } else {
00309 throw error_msg("DLL_ERROR", "Unknown DLL type");
00310 }
00311 (*this).change_ind = clock();
00312
00313 }
00314
00317 void DiffusionCoefficient::MakeDLL_B(GridElement &L, GridElement &pc, GridElement &alpha, double Kp) {
00318 int il, im, ia;
00319 for (il = 0; il < L.size; il++) {
00320 for (im = 0; im < pc.size; im++) {
00321 for (ia = 0; ia < alpha.size; ia++) {
00322 (*this)[il][im][ia] = VF::Df(L[il][im][ia], Kp);
00323 }
00324 }
00325 }
00326 }
00327
00330 void DiffusionCoefficient::MakeDLL_FAKE(GridElement &L, GridElement &pc, GridElement &alpha, double Kp) {
00331 int il, im, ia;
00332 for (il = 0; il < L.size; il++) {
00333 for (im = 0; im < pc.size; im++) {
00334 for (ia = 0; ia < alpha.size; ia++) {
00335 (*this)[il][im][ia] = VF::Df(L[il][im][ia], Kp) + 3e3*pow(10.0, 0.506*Kp-9.325)*pow(L[il][im][ia],5);;
00336 }
00337 }
00338 }
00339 }
00340
00343 void DiffusionCoefficient::MakeDLL_BE(GridElement &L, GridElement &pc, GridElement &alpha, double Kp) {
00344 int il, im, ia;
00345 double e, T, E0, M, E, wd, Re;
00346
00347 e = 4.8e-10 * 3.33564e-10;
00348 Re = 6.4e8;
00349 T = 0.75;
00350 E0 = 0.511;
00351
00352 for (il = 0; il < L.size; il++) {
00353 for (im = 0; im < pc.size; im++) {
00354 for (ia = 0; ia < alpha.size; ia++) {
00355 M = VF::pc2mu(L[il][im][ia], pc[il][im][ia], alpha[il][im][ia]);
00356
00357 E = 0.26 * (Kp-1) + 0.1;
00358
00359 wd = ( 3.0 * M * VC::c / (e / pow(L[il][im][ia]*Re, 2) ) ) / sqrt(1 + 2*M*VF::B(L[il][im][ia]) / E0);
00360 (*this)[il][im][ia] = 0.25 * pow(( VC::c * E / VF::B(L[il][im][ia]) ), 2) * ( T / (1 + pow((wd * T / 2), 2)) ) * pow(L[il][im][ia], 6);
00361
00362 (*this)[il][im][ia] = (*this)[il][im][ia]/2;
00363
00364
00365 (*this)[il][im][ia] = (*this)[il][im][ia] + VF::Df(L[il][im][ia], Kp);
00366 }
00367 }
00368 }
00369 }
00370
00371
00374 void DiffusionCoefficient::MakeDLL100(GridElement &L, GridElement &pc, GridElement &alpha, double Kp) {
00375 int il, im, ia;
00376
00377 for (il = 0; il < L.size; il++) {
00378 for (im = 0; im < pc.size; im++) {
00379 for (ia = 0; ia < alpha.size; ia++) {
00380 (*this)[il][im][ia] = VF::max(VF::Df(L[il][im][ia], Kp), 1.0/100);
00381 }
00382 }
00383 }
00384
00385 }
00386
00389 DiffusionCoefficient& DiffusionCoefficient::operator= (const Matrix3D<double> &M) {
00390 Matrix3D<double>::operator=(M);
00391 return *this;
00392 }
00393
00396 DiffusionCoefficient& DiffusionCoefficient::operator= (double val) {
00397 Matrix3D<double>::operator=(val);
00398 return *this;
00399 }
00400
00403 DiffusionCoefficient DiffusionCoefficient::operator* (double val) {
00404 DiffusionCoefficient tmp(*this);
00405 tmp = tmp.Matrix3D<double>::operator*(val);
00406 return tmp;
00407 }
00408
00411 DiffusionCoefficient DiffusionCoefficient::operator+ (Matrix3D<double> &M) {
00412 DiffusionCoefficient tmp(*this);
00413 tmp = tmp.Matrix3D<double>::operator+(M);
00414 return tmp;
00415 }
00416
00419 DiffusionCoefficient DiffusionCoefficient::operator* (Matrix3D<double> &M) {
00420 DiffusionCoefficient tmp(*this);
00421 tmp = tmp.Matrix3D<double>::operator*(M);
00422 return tmp;
00423 }
00424
00425
00426
00430 DiffusionCoefficientsGroup::DiffusionCoefficientsGroup(string type, Grid &grid) : DiffusionCoefficient(grid) {
00431
00432 this->type = type;
00433
00434 initialized = false;
00435 }
00436
00439 DiffusionCoefficientsGroup::DiffusionCoefficientsGroup(string type, int L_size, int pc_size, int alpha_size) : DiffusionCoefficient(L_size, pc_size, alpha_size) {
00440
00441 this->type = type;
00442
00443 initialized = false;
00444 }
00445
00448 DiffusionCoefficientsGroup::DiffusionCoefficientsGroup(string type, Grid &grid, DiffusionCoefficientParamStructureList DxxParamStructureList) : DiffusionCoefficient(grid){
00449
00450 this->type = type;
00451
00452 Get(grid, DxxParamStructureList);
00453 }
00454
00455
00456
00461 void DiffusionCoefficientsGroup::Get(Grid &grid, DiffusionCoefficientParamStructureList DxxParamStructureList) {
00462
00463 Get(grid, DxxParamStructureList, this->type);
00464 }
00465
00472 void DiffusionCoefficientsGroup::Get(Grid &grid, DiffusionCoefficientParamStructureList DxxParamStructureList, string type) {
00473 unsigned int Dxx_it;
00474 for (Dxx_it = 0; Dxx_it < DxxParamStructureList.size(); Dxx_it++) {
00475
00476 if (DxxParamStructureList[Dxx_it].DxxType == type) {
00477
00478 this->name = DxxParamStructureList[Dxx_it].DxxName;
00479
00480 DxxList.push_back(DiffusionCoefficient(grid, DxxParamStructureList[Dxx_it]));
00481 this->change_ind = clock();
00482 } else {
00483
00484 }
00485 }
00486
00487 }
00488
00489
00493 bool DiffusionCoefficientsGroup::ActivateAndScale(double time, double Kp) {
00494
00495
00496 double scaleCoeff = 0;
00497 unsigned int Dxx_it;
00498
00499
00500 bool need_changes = false;
00501 for (Dxx_it = 0; Dxx_it < DxxList.size(); Dxx_it++) {
00502 if (time >= DxxList[Dxx_it].Dxx_parameters.time_start && time < DxxList[Dxx_it].Dxx_parameters.time_end) {
00503
00504 if (!DxxList[Dxx_it].is_active) {
00505 need_changes = true;
00506 break;
00507 }
00508
00509 if (DxxList[Dxx_it].useScale) {
00510 need_changes = true;
00511 break;
00512 }
00513 } else {
00514
00515 if (DxxList[Dxx_it].is_active) {
00516 need_changes = true;
00517 break;
00518 }
00519 }
00520 }
00521 if (need_changes == false) return false;
00522
00523 (*this) = 0;
00524 this->change_ind = clock();
00525
00526 for (Dxx_it = 0; Dxx_it < DxxList.size(); Dxx_it++) {
00527
00528 if (time >= DxxList[Dxx_it].Dxx_parameters.time_start && time < DxxList[Dxx_it].Dxx_parameters.time_end) {
00529
00530 if (!DxxList[Dxx_it].is_active) {
00531 DxxList[Dxx_it].is_active = true;
00532 Output::echo(2, "Activeted %s.\n", DxxList[Dxx_it].name.c_str());
00533 }
00534
00535 if (DxxList[Dxx_it].useScale) {
00536
00537 scaleCoeff = DxxList[Dxx_it].Scale(Kp);
00538 Output::echo(2, "Bw scaled %s by %f.\n", DxxList[Dxx_it].name.c_str(), scaleCoeff);
00539
00540 (*this) = DxxList[Dxx_it] * scaleCoeff + (*this);
00541 } else {
00542 (*this) = (*this) + DxxList[Dxx_it];
00543 }
00544 } else {
00545
00546 if (DxxList[Dxx_it].is_active) {
00547 DxxList[Dxx_it].is_active = false;
00548 Output::echo(2, "Deactivated %s.", DxxList[Dxx_it].name.c_str());
00549 }
00550 }
00551 }
00552 return false;
00553 }
00554
00558 double DiffusionCoefficient::Scale(double Kp) {
00559 double Bw2=0, Bw2_0=0;
00560 if (!this->useScale) return 1;
00561 if (Dxx_parameters.waveType == "WT_CHORUS_DAY" || Dxx_parameters.waveType == "WT_CHORUS_NIGHT") {
00562
00563 if (Dxx_parameters.DxxKp <= 2.333) Bw2_0 = 2.0*pow(10, 0.73 + 0.91 * Dxx_parameters.DxxKp);
00564 if (Dxx_parameters.DxxKp > 2.333) Bw2_0 = 2.0*pow(10, 2.5 + 0.18 * Dxx_parameters.DxxKp);
00565 if (Kp <= 2.333) Bw2 = 2.0*pow(10, 0.73 + 0.91 * Kp);
00566 if (Kp > 2.333) Bw2 = 2.0*pow(10, 2.5 + 0.18 * Kp);
00567 return (Bw2/Bw2_0);
00568 } else if (Dxx_parameters.waveType == "WT_HISS_PP" || Dxx_parameters.waveType == "WT_HISS_PLUME") {
00569
00570 return pow(Kp/Dxx_parameters.DxxKp, 2);
00571 } else if (Dxx_parameters.waveType == "WT_EMIC_PLUME" ) {
00572
00573 return pow(Kp/Dxx_parameters.DxxKp, 4);
00574 } else {
00575
00576 return 1;
00577 }
00578 }
00579
00583 DiffusionCoefficientsGroup& DiffusionCoefficientsGroup::operator= (double val) {
00584
00585 DiffusionCoefficient::operator=(val);
00586 return *this;
00587 }
00588
00592 DiffusionCoefficientsGroup& DiffusionCoefficientsGroup::operator= (const Matrix3D<double> &M) {
00593 DiffusionCoefficient::operator=(M);
00594 return *this;
00595 }
00596
00600 DiffusionCoefficientsGroup& DiffusionCoefficientsGroup::operator= (const DiffusionCoefficient &Dxx) {
00601 DiffusionCoefficient::operator=(Dxx);
00602 return *this;
00603 }
00604
00605
00606
00607
00608
00609 double Alpha_ne(double pangle, double lambda, double L) {
00610 return asin( sqrt( f1(lambda)) * sin(pangle) );
00611 }
00612
00613 double f1(double lambda) {
00614 return sqrt(4.0 - 3.0 * pow(cos(lambda), 2))/(pow(cos(lambda), 6));
00615 }
00616
00617 double B(double lambda, double L) {
00618 return f1(lambda) / pow(L, 3) * VC::B_0;
00619 }
00620
00621 double func_tmp (double x, double Alpha) {
00622 return pow(x, 6) + (3.0 * pow(sin(Alpha), 4)) * x - 4.0 * pow(sin(Alpha), 4);
00623 }
00624
00625 double F_cap(double x, double y, double b, double s, double epsilon, DiffusionCoefficientParamStructure DxxParamStructure) {
00626 double c1 = 2.0*s*(-1+epsilon);
00627 double c2 = 1.0-4.0*epsilon+pow(epsilon,2);
00628 double c3 = -s*(-1+epsilon)*(b+4.0*epsilon)/2.0;
00629 double c4 = epsilon*(b+epsilon) ;
00630 double g = pow(x,4) + c1*pow(x,3) + c2*pow(x,2) + c3*x + c4;
00631 return y*(pow((x-s),2))*(pow((x+s*epsilon),2))/(x*g);
00632 }
00633
00634 double F_cap2(double x, double y, double a, double beta, double mu, double s, double epsilon, double Alpha_star, DiffusionCoefficientParamStructure DxxParamStructure) {
00635 double up = 2.0*(x+a)/(beta*mu);
00636 double down = 2.0*x-1.0/Alpha_star*(1.0/(x-s)+DxxParamStructure.eta1*epsilon/(x+s*epsilon)+DxxParamStructure.eta2*epsilon/(4.0*x+s*epsilon)+DxxParamStructure.eta3*epsilon/(16.0*x+s*epsilon))+
00637 x*1.0/Alpha_star*(1.0/(pow((x-s), 2))+DxxParamStructure.eta1*epsilon/(pow((x+s*epsilon), 2))+4.0*epsilon*DxxParamStructure.eta2/(pow((4.0*x+s*epsilon), 2))+16.0*epsilon*DxxParamStructure.eta3/(pow((16.0*x+s*epsilon), 2)));
00638 return up/down;
00639 }
00640
00641 double quad1(double (*func)(double lambda, DiffusionCoefficientParamStructure DxxParamStructure), double a, double b, int M, DiffusionCoefficientParamStructure DxxParamStructure) {
00642
00643
00644
00645 double h = (b - a)/M;
00646 double s = func(a, DxxParamStructure);
00647 double x;
00648 int k;
00649 for (k = 1; k <= (M-1); k++) {
00650 x = a + h*k;
00651 s = s + func(x, DxxParamStructure);
00652 }
00653 s = h * (func(a, DxxParamStructure) + func(b, DxxParamStructure))/2.0 + h*s;
00654 return s;
00655 }
00656
00657 double Dxx_ba(double L1, double EMeV, double Alpha, double int_Dxx_loc(double lambda, DiffusionCoefficientParamStructure DxxParamStructure), DiffusionCoefficientParamStructure DxxParamStructure) {
00658
00659
00660 DxxParamStructure.L = L1;
00661 DxxParamStructure.EMeV = EMeV;
00662 DxxParamStructure.Alpha = Alpha;
00663
00664 double lambda_m;
00665
00666 lambda_m = acos(sqrt(bisection (func_tmp, DxxParamStructure.Alpha, 0, 1)))*0.99;
00667 if (lambda_m > DxxParamStructure.lam_max*VC::pi/180) lambda_m = DxxParamStructure.lam_max*VC::pi/180.;
00668
00669 return (1.0/(1.38-0.32*(sin(DxxParamStructure.Alpha) + pow(sin(DxxParamStructure.Alpha),0.5))) * quad1(int_Dxx_loc, DxxParamStructure.lam_min*VC::pi/180., lambda_m, DxxParamStructure.nint, DxxParamStructure));
00670 }
00671
00672
00673 double int_Daa_loc(double lambda, DiffusionCoefficientParamStructure DxxParamStructure) {
00674
00675 double Daa = Dxx_local(lambda, Daa_root, DxxParamStructure);
00676 if (DxxParamStructure.Alpha*180/VC::pi < 0.001) Daa = 0;
00677 return Daa * cos(Alpha_ne(DxxParamStructure.Alpha, lambda, DxxParamStructure.L)) * pow(cos(lambda), 7) /(pow(cos(DxxParamStructure.Alpha), 2));
00678 }
00679
00680 double int_Dpp_loc(double lambda, DiffusionCoefficientParamStructure DxxParamStructure) {
00681
00682 double Dpp = Dxx_local(lambda, Dpp_root, DxxParamStructure);
00683 if (DxxParamStructure.Alpha*180/VC::pi < 0.001) Dpp = 0;
00684 return Dpp * (cos(lambda) * pow((1 + 3*pow(sin(lambda), 2)), 0.5)) / cos(Alpha_ne(DxxParamStructure.Alpha, lambda, DxxParamStructure.L));
00685 }
00686
00687 double int_Dpa_loc(double lambda, DiffusionCoefficientParamStructure DxxParamStructure) {
00688
00689 double Dpa = Dxx_local(lambda, Dpa_root, DxxParamStructure);
00690 if (DxxParamStructure.Alpha*180/VC::pi < 0.001) Dpa = 0;
00691 return Dpa * sin(Alpha_ne(DxxParamStructure.Alpha, lambda, DxxParamStructure.L)) * pow(cos(lambda), 7) / (cos(DxxParamStructure.Alpha) * sin(DxxParamStructure.Alpha));
00692 }
00693
00694 double Dxx_local(double lambda, double Dxx_root(double Omega_e, double x, double mu, double su, double y, double beta, double a, double b, double Alpha_star, double s, double epsilon, double d_x, double x_m, double R, DiffusionCoefficientParamStructure DxxParamStructure), DiffusionCoefficientParamStructure DxxParamStructure ) {
00695
00696 double Omega_e = VC::qe * B(lambda, DxxParamStructure.L) / VC::m / VC::c;
00697 double Omega_e_eq = VC::qe * B(0, DxxParamStructure.L) / VC::m / VC::c;
00698
00699
00700
00701 double Omega_p;
00702
00703 double N_0;
00704 double LT;
00705 if (DxxParamStructure.numberDensity == "ND_CONSTANT") {
00706 Omega_p = DxxParamStructure.f * Omega_e_eq;
00707 } else if (DxxParamStructure.numberDensity == "ND_CHORUS") {
00708 N_0 = 124.0*(pow((3.0/DxxParamStructure.L),4));
00709 Omega_p = sqrt(4.*VC::pi*N_0*pow(VC::qe,2)/VC::m);
00710 DxxParamStructure.f = Omega_p / Omega_e_eq;
00711
00712 } else if (DxxParamStructure.numberDensity == "ND_CHORUS_NIGHT") {
00713 LT = 3;
00714 N_0 = 124.0*(pow((3.0/DxxParamStructure.L),4)) + 36.0*pow((3.0/DxxParamStructure.L),3.5) * cos((LT-(7.7*pow((3.0/DxxParamStructure.L),2.0)+12))*VC::pi/12);
00715 Omega_p = sqrt(4.*VC::pi*N_0*pow(VC::qe,2)/VC::m);
00716 DxxParamStructure.f = Omega_p / Omega_e_eq;
00717
00718 } else if (DxxParamStructure.numberDensity == "ND_CHORUS_DAY") {
00719 LT = 9;
00720 N_0 = 124.0*(pow((3.0/DxxParamStructure.L),4)) + 36.0*pow((3.0/DxxParamStructure.L),3.5) * cos((LT-(7.7*pow((3.0/DxxParamStructure.L),2.0)+12))*VC::pi/12);
00721 Omega_p = sqrt(4.*VC::pi*N_0*pow(VC::qe,2)/VC::m);
00722 DxxParamStructure.f = Omega_p / Omega_e_eq;
00723
00724 } else if (DxxParamStructure.numberDensity == "ND_HISS") {
00725 N_0 = pow(10, (-0.3145*DxxParamStructure.L + 3.9043));
00726 Omega_p = sqrt(4.*VC::pi*N_0*pow(VC::qe,2)/VC::m);
00727 DxxParamStructure.f = Omega_p / Omega_e_eq;
00728
00729 } else if (DxxParamStructure.numberDensity == "ND_HISS2") {
00730 N_0 = 1e3*pow(4.0/DxxParamStructure.L, 4);
00731 Omega_p = sqrt(4.*VC::pi*N_0*pow(VC::qe,2)/VC::m);
00732 DxxParamStructure.f = Omega_p / Omega_e_eq;
00733
00734 } else {
00735 throw error_msg("DXX_DENSITY_ERROR", "Unknown number density type: %s\n", DxxParamStructure.numberDensity.c_str());
00736 }
00737
00738
00739 double epsilon = 1.0/1837;
00740 double beta = (pow((DxxParamStructure.EMeV * (DxxParamStructure.EMeV + 2)), 0.5)) / (DxxParamStructure.EMeV + 1);
00741 double mu = cos(Alpha_ne(DxxParamStructure.Alpha, lambda, DxxParamStructure.L));
00742 double su = sin(Alpha_ne(DxxParamStructure.Alpha, lambda, DxxParamStructure.L));
00743
00744 double lam_sp = -1;
00745 double gamma = DxxParamStructure.EMeV + 1;
00746 double a = DxxParamStructure.s * lam_sp / gamma;
00747
00748 double Alpha_star = pow((Omega_e / Omega_p), 2);
00749 double b = (1 + epsilon)/Alpha_star;
00750
00751
00752
00753 double d_omega;
00754 double Omega_m;
00755 double Omega_1;
00756 double Omega_2;
00757
00758 if (DxxParamStructure.Omega_mType == "OT_PARTIAL") {
00759 Omega_m = DxxParamStructure.Omega_m * Omega_e_eq;
00760 d_omega = DxxParamStructure.d_omega * Omega_e_eq ;
00761 Omega_1 = DxxParamStructure.omega_lc * Omega_e_eq;
00762 Omega_2 = DxxParamStructure.omega_uc * Omega_e_eq;
00763 } else if (DxxParamStructure.Omega_mType == "OT_ABSOLUTE") {
00764 Omega_m = DxxParamStructure.Omega_m * 2 * VC::pi;
00765 d_omega = DxxParamStructure.d_omega * 2 * VC::pi;
00766 Omega_1 = DxxParamStructure.omega_lc * 2 * VC::pi;
00767 Omega_2 = DxxParamStructure.omega_uc * 2 * VC::pi;
00768 } else {
00769 throw error_msg("DXX_OMEGA_M_TYPE", "Unknown Omega_m type: ", DxxParamStructure.Omega_mType.c_str());
00770 }
00771
00772
00773
00774
00775
00776
00777
00778 double Bw;
00779 if (DxxParamStructure.BwFromLambda) {
00780 Bw = pow(10, 0.75 + 0.04*(lambda*180/VC::pi))*1.e-8 * DxxParamStructure.Bw;
00781 } else {
00782 Bw = DxxParamStructure.Bw;
00783 }
00784 double R = pow((Bw / B(lambda, DxxParamStructure.L)), 2);
00785
00786 double x_m=0;
00787 double d_x=0;
00788 double x_1=0;
00789 double x_2=0;
00790 if (DxxParamStructure.particle == "PT_ELECTRON") {
00791 x_m = Omega_m / Omega_e;
00792 d_x = d_omega / Omega_e;
00793 x_1 = Omega_1 / Omega_e;
00794 x_2 = Omega_2 / Omega_e;
00795 } else if (DxxParamStructure.particle == "PT_PROTON") {
00796 double cfm = 16;
00797 x_m = Omega_m / Omega_e * epsilon / cfm;
00798 d_x = d_omega / Omega_e * epsilon / cfm;
00799 x_1 = Omega_1 / Omega_e * epsilon / cfm;
00800 x_2 = Omega_2 / Omega_e * epsilon / cfm;
00801 }
00802
00803 std::vector<double> xx = rrouts(x_1, x_2, DxxParamStructure.eta1, DxxParamStructure.eta2, DxxParamStructure.eta3, epsilon, beta, mu, Alpha_star, a, DxxParamStructure);
00804
00805 double x, y;
00806 double Dxx = 0, Dxx_tmp = 0;
00807 for (unsigned int i = 1; i <= xx.size();i++) {
00808 x = xx[i-1];
00809 y = (x + a) / (beta * mu);
00810 Dxx_tmp = Dxx_root(Omega_e, x, mu, su, y, beta, a, b, Alpha_star, DxxParamStructure.s, epsilon, d_x, x_m, R, DxxParamStructure);
00811 Dxx = Dxx_tmp + Dxx;
00812 }
00813 if (xx.size() == 0) {
00814 Dxx = 0;
00815 }
00816
00817 return Dxx;
00818 }
00819
00820 double Daa_root(double Omega_e, double x, double mu, double su, double y, double beta, double a, double b, double Alpha_star, double s, double epsilon, double d_x, double x_m, double R, DiffusionCoefficientParamStructure DxxParamStructure) {
00821 return (VC::pi/2) * (1.0/DxxParamStructure.nu) * Omega_e * (1.0/(pow((DxxParamStructure.EMeV + 1), 2))) *
00822 ((R * pow((1 - x * mu / (y * beta)), 2) * fabs(F_cap2(x,y,a,beta,mu,s,epsilon,Alpha_star, DxxParamStructure)))/
00823 (d_x * fabs(beta * mu - F_cap2(x,y,a,beta,mu,s,epsilon,Alpha_star, DxxParamStructure))))*
00824 exp( - (pow(((x - x_m) / d_x), 2)));
00825 }
00826
00827 double Dpa_root(double Omega_e, double x, double mu, double su, double y, double beta, double a, double b, double Alpha_star, double s, double epsilon, double d_x, double x_m, double R, DiffusionCoefficientParamStructure DxxParamStructure) {
00828 return -(VC::pi/2)*(1.0/DxxParamStructure.nu)* Omega_e*su/beta*(1.0/(pow((DxxParamStructure.EMeV+1), 2)))*
00829 ((R*(1-x/y*mu/beta)*(x/y)* fabs(F_cap2(x,y,a,beta,mu,s,epsilon,Alpha_star,DxxParamStructure)))/
00830 (d_x * fabs(beta * mu - F_cap2(x,y,a,beta,mu,s,epsilon,Alpha_star, DxxParamStructure))))*
00831 exp( -(pow(((x - x_m) / d_x), 2)));
00832 }
00833
00834 double Dpp_root(double Omega_e, double x, double mu, double su, double y, double beta, double a, double b, double Alpha_star, double s, double epsilon, double d_x, double x_m, double R, DiffusionCoefficientParamStructure DxxParamStructure) {
00835 return (VC::pi/2) * (1.0/DxxParamStructure.nu) * (1./pow(beta, 2)) * (1-pow(mu, 2)) * Omega_e * (1.0/(pow((DxxParamStructure.EMeV + 1), 2))) *
00836
00837
00838 ((R * (pow(x/y, 2)) * fabs(F_cap2(x,y,a,beta,mu,s,epsilon,Alpha_star,DxxParamStructure)))/
00839 (d_x*fabs(beta *mu-F_cap2(x,y,a,beta,mu,s,epsilon,Alpha_star, DxxParamStructure))))*
00840 exp(-(pow(((x-x_m)/d_x),2)));
00841 }
00842
00843
00844
00845 std::vector<double> rrouts(double x_1, double x_2, double eta1, double eta2, double eta3, double epsilon, double beta, double mu, double Alpha_star, double a, DiffusionCoefficientParamStructure DxxParamStructure) {
00846
00847
00848
00849
00850
00851
00852 double c0=-DxxParamStructure.s;
00853 double c1=DxxParamStructure.s*epsilon;
00854 double c2=1.0/4.*DxxParamStructure.s*epsilon;
00855 double c3=1.0/16.*DxxParamStructure.s*epsilon;
00856
00857 double d0=1;
00858 double d1=eta1*epsilon;
00859 double d2=eta2/4.*epsilon;
00860 double d3=eta3/16.*epsilon;
00861
00862 double g3=d0+d1+d2+d3;
00863 double g2=d0*(c1+c2+c3)+d1*(c0+c2+c3)+d2*(c0+c1+c3)+d3*(c0+c1+c2);
00864 double g1=d0*(c1*c2+c1*c3+c2*c3)+d1*(c0*c2+c0*c3+c2*c3)+d2*(c0*c1+c0*c3+c1*c3)+d3*(c0*c1+c0*c2+c1*c2);
00865 double g0=d0*c1*c2*c3+d1*c0*c2*c3+d2*c0*c1*c3+d3*c0*c1*c2;
00866
00867 double h4=1.;
00868 double h3=c0+c1+c2+c3;
00869 double h2=c0*c1+c0*c2+c0*c3+c1*c2+c1*c3+c2*c3;
00870 double h1=c0*c1*c2+c0*c1*c3+c0*c2*c3+c1*c2*c3;
00871 double h0=c0*c1*c2*c3;
00872
00873 double A[7];
00874 A[0] = (pow(beta,2)*pow(mu,2)-1.)*h4;
00875 A[1] = (pow(beta,2)*pow(mu,2)-1.)*h3-2*a*h4;
00876 A[2] = (pow(beta,2)*pow(mu,2)-1.)*h2-2*a*h3-pow(beta,2)*pow(mu,2)*1./Alpha_star*g3-pow(a,2)*h4;
00877 A[3] = (pow(beta,2)*pow(mu,2)-1.)*h1-2*a*h2-pow(beta,2)*pow(mu,2)*1./Alpha_star*g2-pow(a,2)*h3;
00878 A[4] = (pow(beta,2)*pow(mu,2)-1.)*h0-2*a*h1-pow(beta,2)*pow(mu,2)*1./Alpha_star*g1-pow(a,2)*h2;
00879 A[5] = -2*a*h0-pow(beta,2)*pow(mu,2)*1./Alpha_star*g0-pow(a,2)*h1;
00880 A[6] = -pow(a,2)*h0;
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890 std::vector<double> xx;
00891 double r[7], wr[6],wi[6],quad[2];
00892 int n, numr;
00893 n=6;
00894
00895 quad[0] = 2.71828e-1;
00896 quad[1] = 3.14159e-1;
00897
00898 get_quads(A, n, quad, r);
00899 numr = roots(r, n, wr, wi);
00900
00901 for (int i = 0; i < numr; i++) {
00902 if (fabs(wi[i]) < double_zero && wr[i] < x_2 && wr[i] > x_1) {
00903 xx.push_back(wr[i]);
00904 }
00905 }
00906
00907 return xx;
00908
00909 }