/*************************************************************************** * $Id: models21to30.cpp,v 1.8 2009/03/17 14:44:43 thernis Exp $ * ****************************************************************************/ #include "models21to30.h" #include #include #include "Cvec.h" #include "constant.h" #include "rtmiscfunc.h" // -- density 21 // Sphericaly symmetric Ne // // Use the Saito equatorial density model float CModel21::Density(const Cvec &v) { float r=v.norm(); if (r <= 1.05) return 0.; // -- coeff for the equatorial density model float c1=1.36e6; float d1=2.14; float c2=1.68e8; float d2=6.13; float nel=c1*pow(r,-d1)+c2*pow(r,-d2); return nel; } void CModel21::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0x4; vp.push_back(moddefparam("","Saito equatorial density.","","")); return; } // -- density 22 // Same as density 14 but with a modified value of the angular sector // see 11 May 2004 float CModel22::Density(const Cvec &v) { float x=v.v[0],y=v.v[1],z=v.v[2],r=v.norm(); if (r <= 2.57) return 0.; float rbidon,phi,theta; cvcoord(x,y,z,&rbidon,&phi,&theta); //if (phi < 1.0472 or phi > 2.0944) return 0.; // if (phi < (45*DTOR) or phi > (135*DTOR)) return 0.; if (phi < (45*DTOR) or phi > (135*DTOR)) 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 // see 7 May 2004 // *2 correction fudge factor (don't know why that works !) // -- 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 cw[4]={17.79/4.71,-85.93/4.71,138.32/4.71,391.45/4.71}; // valid from 2.5 to 20 Rsun float cw[4]={17.79*DTOR*DTOR,-85.93*DTOR*DTOR,138.32*DTOR*DTOR,391.45*DTOR*DTOR}; // valid from 2.5 to 20 Rsun // see 10 May 2004 float cu[4]={60.9441*DTOR , -542.687*DTOR , 1889.58*DTOR , -2289.55*DTOR}; float coef[4]={-1,-2,-3,-4}; float nel=0; // -- streamer half thickness float w0=0,u0=0; for(int i=0;i<=3;i++) { float powr=pow(r,coef[i]); nel+=c[i]*powr; w0+=cw[i]*powr; u0+=cu[i]*powr; } float dcross=w0/(2*u0); if (fabs(theta) < dcross) { float ee=(theta*theta)/w0; if (ee > 1E1) return 0.; nel*=exp(-ee); } else { float ee=fabs(theta/u0)-dcross/u0/2; if (ee > 1E1) return 0.; nel*=exp(-ee); } return nel; } void CModel22::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0x4; vp.push_back(moddefparam("","Streamer slab model, with a modified value of the angular sector compare to model 14.","","")); return; } // -- density 23 // Same as density 18 but with a modified value of the angular sector // see 11 May 2004 // see 14 June 2004 float CModel23::Density(const Cvec &v) { float x=v.v[0],y=v.v[1],z=v.v[2],r=v.norm(); if (r <= 2.57) return 0.; float rbidon,phi,theta; cvcoord(x,y,z,&rbidon,&phi,&theta); float sechang=*(pparam+4); if (phi < (PI/2-sechang) or phi > (PI/2+sechang)) 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 // see 7 May 2004 // *2 correction fudge factor (don't know why that works !) // float cw[4]={17.79/4.71,-85.93/4.71,138.32/4.71,391.45/4.71}; // valid from 2.5 to 20 Rsun float cw[4]={17.79*DTOR*DTOR,-85.93*DTOR*DTOR,138.32*DTOR*DTOR,391.45*DTOR*DTOR}; // valid from 2.5 to 20 Rsun // see 10 May 2004 float cu[4]={60.9441*DTOR , -542.687*DTOR , 1889.58*DTOR , -2289.55*DTOR}; float coef[4]={-1,-2,-3,-4}; float nel=0; // -- streamer half thickness float w0=0,u0=0; float *pp; pp=pparam; for(int i=0;i<=3;i++) { float powr=pow(r,coef[i]); // nel+=c[i]*powr; nel+=(*(pp++))*powr; w0+=cw[i]*powr; u0+=cu[i]*powr; } float dcross=w0/(2*u0); if (fabs(theta) < dcross) { float ee=(theta*theta)/w0; if (ee > 1E1) return 0.; nel*=exp(-ee); } else { float ee=fabs(theta/u0)-dcross/u0/2; if (ee > 1E1) return 0.; nel*=exp(-ee); } return nel; } void CModel23::initParam(float *pparam) { this->pparam=pparam; } void CModel23::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0; vp.push_back(moddefparam("","Streamer slab model, with a modified value of the angular sector compare to model 18.","","")); vp.push_back(moddefparam("c","[1.34e4,1.15e6,-6.022e6,5.577e7,0.698]","Coefficients for the Hayes model.","")); return; } // -- density 24 // Cylinder model, on axis Ox float CModel24::Density(const Cvec &v) { float x=v.v[0],y=v.v[1],z=v.v[2],r=v.norm(); if (r <= 2.5 || x <= 1.1) return 0.; float rbidon,phi,theta; cvcoord(x,y,z,&rbidon,&phi,&theta); float dist=sqrt(z*z+y*y); //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.2; for(int i=0;i<=3;i++) { nel+=c[i]*pow(r,coef[i]); } //nel=1./r; float ee=dist/d0; ee*=ee; if (ee > 1E1) return 0.; nel*=exp(-ee); return nel; } void CModel24::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0x4; vp.push_back(moddefparam("","Cylinder model, Hayes radial, exp(-abs(d/d0)) orthoradial, axis Ox","","")); return; } // -- density 25 // Datacube passed in pparam pointer float CModel25::Density(const Cvec &v) { // -- nothing inside the Sun ! if (v.norm() <= 1.01) return 0.; // ---- compute position of the voxel in the data cube // return 0 if out of the cube // -- X float xf=(pparam[3]+v.v[0]/voxsize); if (xf < 0 || xf > (pparam[0])) return 0; long xp=long(floor(double(xf))); // -- Y float yf=(pparam[4]+v.v[1]/voxsize); if (yf < 0 || yf > (pparam[1])) return 0; long yp=long(floor(double(yf))); // -- Z float zf=(pparam[5]+v.v[2]/voxsize); if (zf < 0 || zf > (pparam[2])) return 0; long zp=long(floor(double(zf))); // -- position of the voxel: long offset=7 + xp + yp*sx + zp*sx*sy; return pparam[offset]; } void CModel25::initParam(float *pparam) { this->pparam=pparam; // -- get the voxel side size in rsun voxsize=pparam[6]; sx=long(pparam[0]); sy=long(pparam[1]); } void CModel25::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0; vp.push_back(moddefparam("","Density cube.","","")); vp.push_back(moddefparam("cubefile","getenv('RT_PATH')+get_delim()+'testcube.fts'","Filename of the density cube","")); vp.push_back(moddefparam("mp","","","")); vp.push_back(moddefparam("","mp=loaddenscube($cubefile)","Load the density cube.","")); return; } // -- density 26 // Datacube passed in pparam pointer. float CModel26::Density(const Cvec &vv) { // -- compute position of the voxel in the data cube // return 0 if out of the cube // ---- X float xf=( (vv[0]/voxsize) + pparam[3]); if (xf < 0 || xf >=pparam[0]) return 0; float xs,t; long xp; if (xf <= 0.5) { xs=0.;xp=0;t=0.; } else if (xf >= (pparam[0]-0.5)) { xs=pparam[0]-2;xp=long(xs);t=1.; } else { xs=xf-0.5;xp=long(xs);t=(xs-float(xp)); } // ---- Y float yf=((vv[1]/voxsize) + pparam[4]); if (yf < 0 || yf >=pparam[1]) return 0; float ys,u; long yp; if (yf <= 0.5) { ys=0.;yp=0;u=0.; } else if (yf >= (pparam[1]-0.5)) { ys=pparam[1]-2;yp=long(ys);u=1.; } else { ys=yf-0.5;yp=long(ys);u=(ys-float(yp)); } // ---- Z float zf=((vv[2]/voxsize) + pparam[5]); if (zf < 0 || zf >=pparam[2]) return 0; float zs,v; long zp; if (zf <= 0.5) { zs=0.;zp=0;v=0.; } else if (zf >= (pparam[2]-0.5)) { zs=pparam[2]-2;zp=long(zs);v=1.; } else { zs=zf-0.5;zp=long(zs);v=(zs-float(zp)); } float *pcube; pcube=pparam+7; // -- trilinear interpolation float neinterp=trilininterp(t,u,v,xp,yp,zp,pcube,sx,sy); return neinterp; } void CModel26::initParam(float *pparam) { this->pparam=pparam; voxsize=pparam[6]; sx=long(pparam[0]); sy=long(pparam[1]); sxsy=sx*sy; } void CModel26::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0; vp.push_back(moddefparam("","Density cube, trilinear interpolation for smoothing.","","")); vp.push_back(moddefparam("cubefile","getenv('RT_PATH')+get_delim()+'testcube.fts'","Filename of the density cube","")); vp.push_back(moddefparam("mp","","","")); vp.push_back(moddefparam("","mp=loaddenscube($cubefile)","Load the density cube.","")); return; } //------------------------------------------------ // -- density 27 // Spherical Shell Ne density float CModel27::Density(const Cvec &v) { float x=v.v[0],y=v.v[1],z=v.v[2];//,r=v.norm(); float r0 = pparam[0]; // 8; //height of center sphere float d = pparam[1]; // 2; //radius of sphere float d0 = pparam[2]; //.3; //1/2 thickness of shell float nemin = pparam[3]; // 1e5 //float nemax = 1e6; float rr = sqrt(x*x+y*y+(z-r0)*(z-r0)); //-r0 if (z < r0) return 0.; //float nel=0; if (fabs(rr-d) <= d0){ return nemin;// + (nemax-nemin)*(z-r0-d-d0)/(2*d+2*d0); } return 0. ; } void CModel27::initParam(float *pparam) { this->pparam=pparam; } void CModel27::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0; vp.push_back(moddefparam("","Spherical Shell.","","")); vp.push_back(moddefparam("r0","8.","height of center sphere","Rsun")); vp.push_back(moddefparam("d","2.","radius of sphere","Rsun")); vp.push_back(moddefparam("d0","0.3","1/2 thickness of shell","Rsun")); vp.push_back(moddefparam("nemin","1e5","Ne","cm^-3")); return; } //----------------------------------------------- // -- density 28 // Cylindrical Shell Ne density float CModel28::Density(const Cvec &v) { float x=v.v[0],y=v.v[1],z=v.v[2];//,r=v.norm(); float r0 = pparam[0]; //8; //height of center sphere float d = pparam[1]; //2; //radius of sphere float len = pparam[2]; //2; //length of cylender in each direction float d0 = pparam[3]; //0.3; //1/2 thickness of shell float nemin = 1; float nemax = 1; float rr = sqrt(x*x+(z-r0)*(z-r0)); float nel=0; if (fabs(y) <= len) if (fabs(rr-d) <= d0) { nel = nemin + (nemax-nemin) * (z-r0-d-d0)/(2*d+2*d0); //ofstream myFile("output22.txt",ios::app); //myFile<pparam=pparam; } void CModel28::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0; vp.push_back(moddefparam("","Cylindrical Shell.","","")); vp.push_back(moddefparam("r0","8.","height of axis","Rsun")); vp.push_back(moddefparam("d","2.","radius of shell","Rsun")); vp.push_back(moddefparam("len","2.","Lenght of the cylinder","Rsun")); vp.push_back(moddefparam("d0","0.3","Half thickness of the shell","Rsun")); return; } //----------------------------------------------- // -- density 29 // Bow Shock Model by Angelos Vourlidas float CModel29::Density(const Cvec &v) { //pmparam=d,s,h,d0,nel0 float x=v.v[0],y=v.v[1],z=v.v[2];//,r=v.norm(); float d = pparam[0]; //1; float s = pparam[1]; //1.5; float h = pparam[2]; //1;? float d0 = pparam[3]; //.1 float nel0 =pparam[4];//1e5 //s = 1/h; //s/=2; //s+=2; //d+=h/5; float zz = h - (d/s)*pow((sqrt(x*x+y*y)/d),s); float nel = 0; if (fabs(zz-z) <= d0) { nel = nel0; } return nel; } void CModel29::initParam(float *pparam) { this->pparam=pparam; } void CModel29::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0; vp.push_back(moddefparam("","Bow Shock Model (A.Vourlidas)","","")); vp.push_back(moddefparam("d","1.","","Rsun")); vp.push_back(moddefparam("s","1.5","","Rsun")); vp.push_back(moddefparam("h","1","Height of the shell","Rsun")); vp.push_back(moddefparam("d0","0.1","Thickness on the shell","Rsun")); vp.push_back(moddefparam("nel0","1e5","density","cm^-3")); return; } // -- density 30 // Streamer belt simulation with source surface field map: float CModel30::Density(const Cvec &v) { float x=v.v[0],y=v.v[1],z=v.v[2],r=v.norm(); if (r <= 2.55) return 0.; float rbidon,phi,theta; cvcoord(x,y,z,&rbidon,&phi,&theta); // -- param for the neutral sheet map int sang=(int) *pparam,slat=(int) *(pparam+1); // ---- get the distance from the point // to the nearest neighbors // -- nearest neighbor seeking range int srang=15,srlat=15; // in pix // -- point to the neutral sheet map float *pnsheetmap; pnsheetmap=pparam+2; float thetannpos,phinnpos,dist,val; int nnok=wherennsmoothed(pnsheetmap,sang,slat,srang,srlat,(phi*RADEG),(theta*RADEG),&phinnpos,&thetannpos,&dist,&val); if (nnok == 0) return 0.; // -- 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 cw[4]={17.79*DTOR*DTOR,-85.93*DTOR*DTOR,138.32*DTOR*DTOR,391.45*DTOR*DTOR}; // valid from 2.5 to 20 Rsun // see 10 May 2004 float cu[4]={60.9441*DTOR , -542.687*DTOR , 1889.58*DTOR , -2289.55*DTOR}; float coef[4]={-1,-2,-3,-4}; float nel=0; // -- streamer half thickness float w0=0,u0=0; for(int i=0;i<=3;i++) { float powr=pow(r,coef[i]); nel+=c[i]*powr; w0+=cw[i]*powr; u0+=cu[i]*powr; //nel+=c[i]*pow(r,coef[i]); } float dcross=w0/(2*u0); //theta=r*sin(dist*DTOR); // ca c'est tres bourrin ! theta=dist*DTOR; // ca c'est tres bourrin ! //float ee=(theta*theta)/0.005; //float ee=theta/0.01; //if (ee > 1E1) return 0; //nel*=exp(-ee); if (fabs(theta) < dcross) { float ee=(theta*theta)/w0; if (ee > 1E1) return 0.; nel*=val*exp(-ee); } else { float ee=fabs(theta/u0)-dcross/u0/2; if (ee > 1E1) return 0.; nel*=val*exp(-ee); } //nel*=val*exp(-ee); return nel; } void CModel30::initParam(float *pparam) { this->pparam=pparam; } void CModel30::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0; vp.push_back(moddefparam("","Streamer belt simulation with source surface field map, user resolution.","","")); vp.push_back(moddefparam("crot","1912L","Carrington rotation number","")); vp.push_back(moddefparam("modparam","","","")); vp.push_back(moddefparam("","rdtxtmagmap,nsheetmap,crot=$crot","Load PFSS map from WSO web site.","")); vp.push_back(moddefparam("","modparam=reform(360,181,nsheetmap,n_elements(nsheetmap))","Format the array of parameters.","")); return; } /* * $Log: models21to30.cpp,v $ * Revision 1.8 2009/03/17 14:44:43 thernis * Fix bugs in model 26 * * Revision 1.7 2009/02/09 20:51:14 thernis * - Clean up the code * - Change CModel::Density prototype * - Update documentation * - Implement multi-threading using boost thread * - Add new models * * Revision 1.6 2008/09/23 14:08:12 thernis * - add new testing models, - implement integration in front of the instrument * * Revision 1.5 2008/08/18 18:16:15 thernis * Change comments in model 25 * * Revision 1.4 2008/08/18 17:22:11 ontivero * update model 29 * * Revision 1.3 2007/11/20 20:46:32 thernis * Fix rounding problems in model 25 * * Revision 1.2 2007/05/14 17:19:41 thernis * Add CVS id and log in all files * */