/*************************************************************************** * $Id: models01to10.cpp,v 1.3 2009/02/09 20:51:12 thernis Exp $ * ****************************************************************************/ #include "models01to10.h" #include #include "constant.h" #include "rtmiscfunc.h" // -- density 1 // M.Guhathakurta model float CModel01::Density(const Cvec &v,float &temperature) { float x=v.v[0],y=v.v[1],z=v.v[2],r=v.norm(); if (r <= 1.05) return 0.; float alpha=0.*DTOR; float phi0=90.*DTOR; // float nch=1e10; float phi=0,theta; float cp[3]={0.14e7,8.02e7,8.12e7}; float dp[3]={2.8,8.45,16.87}; float ccs[3]={1.41e7,16.42e7,61.90e7}; float dcs[3]={3.39,5.14,12.67}; float gam[3]={16.3,10.0,43.20}; float delt[3]={0.5,7.31,7.52}; float np=0,ncs=0,w=0; float thetamg,nel; theta=asin(z/r); float ratio=x/(r*cos(theta)); if (fabs(ratio) < 1) { phi=acos(ratio); if (y < 0) phi=2.*PI-phi; } for(int i=0;i<=2;i++) { np+=cp[i]*pow(r,-dp[i]); ncs+=ccs[i]*pow(r,-dcs[i]); w+=gam[i]*pow(r,-delt[i]); } thetamg=asin(-cos(theta)*sin(alpha)*sin(phi-phi0)+sin(theta)*cos(alpha)); float expo=pow(thetamg/(w*DTOR),2); if (expo > 1E1) nel=np; else nel=np+(ncs-np)*exp(-expo); return nel; } void CModel01::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0x4; vp.push_back(moddefparam("","M.Guhathakurta model, frozen parameters.","","")); return; } // -- density 2 // Model for testing purpose float CModel02::Density(const Cvec &v,float &temperature) { float x=v.v[0],y=v.v[1],z=v.v[2],r=v.norm(); if (r <= 1.05) return 0.; float nel=0; // float rhosquare=x*x+y*y; if (z > 0.) { float powexp=(x*x+y*y)/(0.5*0.5); if (powexp < 1e1) nel=exp(-r/1.)*exp(-powexp); } return nel; } void CModel02::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0x6; vp.push_back(moddefparam("","Model for testing purpose","","")); return; } // -- density 3 // Model for testing purpose float CModel03::Density(const Cvec &v,float &temperature) { float x=v.v[0],y=v.v[1],z=v.v[2],r=v.norm(); if (r <= 1.05) return 0.; float nel=0; // float theta=0,phi; //phi=acos(z/r); //float ratio=x/(r*sin(phi)); //if (abs(ratio) < 1.) { // theta=acos(ratio); // if (y < 0.) theta=2.*PI-theta; // } float phi,theta; cvcoord(x,y,z,&r,&theta,&phi); //nel=(1./r)*exp(-pow(2.*sin(4*theta),2)/2.)*exp(-pow(phi-PI/2,2)/0.1); if (r > 3 && r < 6 && phi < 0.5 && phi > -0.5 && theta >= 0 && theta <= PI/2) nel=1e4; //nel=(1./r)*exp(-pow(phi-PI/2+0.05*cos(theta),2)/0.01); return nel; } void CModel03::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0x6; vp.push_back(moddefparam("","Model for testing purpose.","","")); return; } // -- density 4 // Density cube from Y.M. with trilinear interpolation float CModel04::Density(const Cvec &v,float &temperature) { return(densyming2(v.v[0],v.v[1],v.v[2],0., pmparam,pmparam+3,pmparam+6,pmparam+9)); } void CModel04::initParam(float* pparam) { this->pmparam=pparam; } void CModel04::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0; vp.push_back(moddefparam("","Density cube from Y-M.Wang, trilinear interpolation in the cube.","","")); vp.push_back(moddefparam("mp","","","")); vp.push_back(moddefparam("","load_dens,dens,rco,phico,thetaco","Load PFSS default density cube.","")); vp.push_back(moddefparam("","mp=[rco,phico,thetaco,reform(dens,n_elements(dens))]","Format the array of parameters.","")); vp.push_back(moddefparam("","dens=0b","Free memory.","")); return; } // -- density 5 //! Same as Model 4 but without trilinear interpolation //! float CModel05::Density(const Cvec &v,float &temperature) { return(densymingnn(v.v[0],v.v[1],v.v[2],0., pmparam,pmparam+3,pmparam+6,pmparam+9)); } void CModel05::initParam(float* pparam) { this->pmparam=pparam; } void CModel05::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0; vp.push_back(moddefparam("","Density cube from Y-M.Wang, no interpolation in the cube.","","")); vp.push_back(moddefparam("modparam","","","")); vp.push_back(moddefparam("","load_dens,dens,rco,phico,thetaco","Load PFSS default density cube.","")); vp.push_back(moddefparam("","modparam=[rco,phico,thetaco,reform(dens,n_elements(dens))]","Format the array of parameters.","")); vp.push_back(moddefparam("","dens=0b","Free memory.","")); return; } // -- density 6 // Very simple model of streamer belt float CModel06::Density(const Cvec &v,float &temperature) { float x=v.v[0],y=v.v[1],z=v.v[2],r=v.norm(); if (r <= 1.05) return 0.; float d0=0.3; // half thickness of the streamer float rbidon,phi,theta; cvcoord(x,y,z,&rbidon,&theta,&phi); float neradial=1/r; float phi0=0.5*cos(2*theta); float d=r*sin(phi-phi0); float neshape=exp(-pow(d/d0,4)); float nel=neradial*neshape; return nel; } void CModel06::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0x6; vp.push_back(moddefparam("","","","")); return; } // -- density 7 // Slab model of streamer // based on Guhathakurta model. // float CModel07::Density(const Cvec &v,float &temperature) { float x=v.v[0],y=v.v[1],z=v.v[2],r=v.norm(); if (r <= 1.05) return 0.; float rbidon,phi,theta; cvcoord(x,y,z,&rbidon,&phi,&theta); if (phi < 0 or phi > (180.*DTOR)) return 0.; // if (abs(theta) > (60./DTOR)) return 0.; // float cp[3]={0.14e7,8.02e7,8.12e7}; // float dp[3]={2.8,8.45,16.87}; float ccs[3]={1.41e7,16.42e7,61.90e7}; float dcs[3]={3.39,5.14,12.67}; float gam[3]={16.3*DTOR,10.0*DTOR,43.20*DTOR}; float delt[3]={0.5,7.31,7.52}; float ncs=0,w=0;//np=0 for(int i=0;i<=2;i++) { //np+=cp[i]*pow(r,-dp[i]); ncs+=ccs[i]*pow(r,-dcs[i]); w+=gam[i]*pow(r,-delt[i]); } // use only the curent sheet model float ee=pow(theta/w,2); if (ee > 1E1) return 0.; float nel=ncs*exp(-ee); //float nel=np; return nel; } void CModel07::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0x6; vp.push_back(moddefparam("","Slab model of streamer, based on Guhathakurta model.","","")); return; } // -- density 8 // Sphericaly symmetric Ne float CModel08::Density(const Cvec &v,float &temperature) { // float x=v.v[0],y=v.v[1],z=v.v[2]; float r=v.norm(); if (r <= 1.05) return 0.; float c[4]={-70544.822 , 1206889.8 , -4264689.2 , 9676682.1}; float coef[4]={-1,-2,-3,-4}; float nel=0; for(int i=0;i<=3;i++) { nel+=c[i]*pow(r,coef[i]); } return nel; } void CModel08::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0x4; vp.push_back(moddefparam("","Uniform corona, Hayes radial.","","")); return; } // -- density 9 // Guhata. radial and Vibert shape float CModel09::Density(const Cvec &v,float &temperature) { float x=v.v[0],y=v.v[1],z=v.v[2],r=v.norm(); if (r <= 1.05) return 0.; float rbidon,phi,theta; cvcoord(x,y,z,&rbidon,&phi,&theta); if (phi < 0 or phi > (180.*DTOR)) return 0.; // if (abs(theta) > (60./DTOR)) return 0.; // float cp[3]={0.14e7,8.02e7,8.12e7}; // float dp[3]={2.8,8.45,16.87}; float ccs[3]={1.41e7,16.42e7,61.90e7}; float dcs[3]={3.39,5.14,12.67}; // float gam[3]={16.3*DTOR,10.0*DTOR,43.20*DTOR}; // float delt[3]={0.5,7.31,7.52}; // -- streamer half thickness float d0=0.2; // float np=0,ncs=0,w=0; float ncs=0; // -- radial term for(int i=0;i<=2;i++) { //np+=cp[i]*pow(r,-dp[i]); ncs+=ccs[i]*pow(r,-dcs[i]); //w+=gam[i]*pow(r,-delt[i]); } // -- shape term //float ee=pow((z/d0),4); float ee=fabs(z/d0); if (ee > 1E1) return 0.; float nel=ncs*exp(-ee); return nel; } void CModel09::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0x6; vp.push_back(moddefparam("","Guhata. radial and Vibert shape","","")); return; } // -- density 10 // Sphericaly symmetric Streamer float CModel10::Density(const Cvec &v,float &temperature) { float x=v.v[0],y=v.v[1],z=v.v[2],r=v.norm(); if (r <= 2.5) return 0.; float rbidon,phi,theta; cvcoord(x,y,z,&rbidon,&phi,&theta); if (phi < 1.0472 or phi > 2.0944) return 0.; //if (y < 0) return 0.; // float c[4]={-1168765.3 , 20505397., -1.0056558e+08 , 1.8379559e+08}; // float c[4]={7.27e3*2,5.899e5*2,-2.781e6*2,2.932e7*2}; // valid from 2.5 to 20 Rsun // -- updated on June 15th using pminimizer02.pro float c[4]={1.34e4,1.15e6,-6.022e6,5.577e7};// valid from 2.5 to 20 Rsun float coef[4]={-1,-2,-3,-4}; //float gam[3]={16.3*DTOR,10.0*DTOR,43.20*DTOR}; //float delt[3]={0.5,7.31,7.52}; //w+=gam[i]*pow(r,-delt[i]); float nel=0; // -- streamer half thickness float d0=0.5; for(int i=0;i<=3;i++) { nel+=c[i]*pow(r,coef[i]); } //nel=1./r; float ee=z/d0; ee*=ee; if (ee > 1E1) return 0.; nel*=exp(-ee); return nel; } void CModel10::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0x6; vp.push_back(moddefparam("","Sphericaly symmetric streamer, Hayes radial, exp(-abs(d/d0)) orthoradial","","")); return; } /* * $Log: models01to10.cpp,v $ * Revision 1.3 2009/02/09 20:51:12 thernis * - Clean up the code * - Change CModel::Density prototype * - Update documentation * - Implement multi-threading using boost thread * - Add new models * * Revision 1.2 2007/05/14 17:19:41 thernis * Add CVS id and log in all files * */