/*************************************************************************** * $Id: models31to40.cpp,v 1.3 2009/02/09 20:51:15 thernis Exp $ * ****************************************************************************/ #include "models31to40.h" #include #include #include #include "Cvec.h" #include "constant.h" #include "rtmiscfunc.h" //----------------------------------------------- // -- density 31 // Graduated Curved Cylindrical Shell Model float CModel31::Density(const Cvec &v) { float x=v.v[0],y=v.v[1],z=v.v[2],r=v.norm(); //COMMON /params/rh,dh,d0,rc,rb,db,nemin,nemax float nel = 0; float nemax = pparam[0]; //1.5; //non-constant density float nemin = pparam[1]; //1; //int makecloud = 0; // Input Parameters: float rh = pparam[2]; //14; // height of leading edge of cylinder at the nose float dh = pparam[3]; //4; // diameter of cylinder at the top float d0 = pparam[4]; //0.3; // 1/2 thickness of shell float rc = pparam[5]; //3; // radius of curvature of the cylinder float rb = pparam[6]; //2; // height of base of cylinder float db = pparam[7]; //.1; // diameter of cylinder at the base float del_d = (dh-db)/(rh-rb); float r0 = rh-dh*0.5-rc; // height of center of curvature of flux rope float r1 = sqrt(r0*r0-(rc+dh*.5)*(rc+dh*.5)); // height of transition of conical legs to curved cylinder float theta = atan2(rc,r1); // half angle of width of curved cylinder // float r2 = r1*cos(theta); // projected height of transition // float temp = rc+d0*.5; // float thetamax = atan2(rc+d0,r1);// max angle including shell // float th = atan2(x,y); // angle of point relative to axis float dely1 = 0; // distance of the point to shell middle on far side float dely2 = 0; // distance of the point to shell middle on near side float delx1 = 0; // same for x float delx2 = 0; // //float rdist = 0; // ;if (y < rb || y > rh) return(nel); // beyond the model if (r < r1) { // in cone region (legs) float rmid = y*tan(theta); // distance of middle of shell to central axis if (x < 0) rmid = -rmid; float rcirc = 0.5*(db+(r-rb)*del_d); // radius of cylinder cross section float rdist = sqrt((x-rmid)*(x-rmid)+z*z); float delr = fabs(rdist-rcirc); if (delr <= d0) { nel=nemin + (r-rb)*(nemax-nemin)/(rh-rb); } // if (makecloud){ // if (nel > 0) { // ofstream myFile("output31.txt",ios::app); // myFile< 0) { // ofstream myFile("output31.txt",ios::app); // myFile< 0) { // ofstream myFile("output31.txt",ios::app); // myFile<pparam=pparam; } void CModel31::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0; vp.push_back(moddefparam("","Graduated Curved Cylindrical Shell Model","","")); vp.push_back(moddefparam("nemax","1.5","Ne max","cm^-3")); vp.push_back(moddefparam("nemin","1.","Ne min","cm^-3")); vp.push_back(moddefparam("rh","14.","Height of leading edge of cylinder at the nose","Rsun")); vp.push_back(moddefparam("dh","4.","Diameter of cylinder at the top","Rsun")); vp.push_back(moddefparam("d0","0.3","Half thickness of shell","Rsun")); vp.push_back(moddefparam("rc","3.","Radius of curvature of the cylinder","Rsun")); vp.push_back(moddefparam("rb","2.","Height of base of cylinder","Rsun")); vp.push_back(moddefparam("db","0.1","Diameter of cylinder at the base","Rsun")); return; } // -- density 32 // Streamer belt simulation with source surface field map: float CModel32::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); // -- streamer thickness parameters // w0=d1, u0=d2 float w0=*(pparam+2),u0=*(pparam+3); // ---- 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+4; 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 is given by the distance of the requested // point to the neutral line theta=dist*DTOR; 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 CModel32::initParam(float *pparam) { this->pparam=pparam; } // -- density 33 // Tube shell model float CModel33::Density(const Cvec &v) { float x=v.v[0],y=v.v[1],z=v.v[2],r=v.norm(); float nel = 0; float d0 = pparam[0]; //0.7; //FULL thickness of shell float rb = pparam[1]; //2.55; //dist to bottom of structure if (r <= rb) return nel; float rbidon,phi,theta; cvcoord(x,y,z,&rbidon,&phi,&theta); // -- Oyz is plane of symmetry x=fabs(x); // -- deal only with first quadrant if (y < 0) return nel; float alpha=pparam[2]*DTOR; //30.*DTOR; // angle between axis and foot [rad] float rf=pparam[3]; //10.; // dist junction line-circle [Rsun] //area between the following two distances from skeleton will have non-zero Ne density //float innerDist = pparam[X]; //1.1 //distance from skeleton to inner radius of shell //float outerDist = pparam[X]; //1.5 //distance from skeleton to outer radius of shell // -- center of the circle depending on the alpha and rf param // rs : dist solar center - shell skeleton circle center float rs=rf*(cos(alpha)+tan(alpha)*sin(alpha)); // rc : radius shell skeleton circle float rc=rf*tan(alpha); float rh = rs+rc; //max height of structure??? //Ne densities at min and max heights (range is from rb to rh) float nemin = 1e4; //nemin=density21(v,pparam,temperature); float nemax = 1e5; nemax=nemin; float ratio = pparam[4]; //ratio of tube radius to height //.2; // -- on which side are we ? foot or shell ? float testside=y+x*tan(alpha)-rs; if (testside < 0) { // -- foot side //return 0.; // Tres bourrin : no feet ! // dist point - foot float r=psldist(Cvec(0,0,0),Cvec(sin(alpha),cos(alpha),0),Cvec(x,y,z)); // constant density for the moment //if (r < d0) return(1e4); //else return(0); if (r < v.norm()*ratio+d0 && r > v.norm()*ratio) { nel = nemin + (v.norm()-rb)*(nemax-nemin)/(rh-rb); } return(nel); //if (r < outerDist && r > innerDist) {return 1e4;} else return 0.; } else { // -- shell side // angular position float rrr,beta,ttt; cvcoord(x,y-rs,0,&rrr,&beta,&ttt); // dist point - circle Cvec vpc(x-rc*cos(beta),y-(rs+rc*sin(beta)),z); // constant density for the moment //if (vpc.norm() < ratio) return(1e4); //else return(0); if (vpc.norm() < v.norm()*ratio+d0 && vpc.norm() > v.norm()*ratio) { nel = nemin + (v.norm()-rb)*(nemax-nemin)/(rh-rb); } return(nel); //if (vpc.norm() < outerDist && vpc.norm() > innerDist) {return 1e4;} else return 0.; } return 0; } void CModel33::initParam(float *pparam) { this->pparam=pparam; } void CModel33::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0; vp.push_back(moddefparam("","Tube shell model.","","")); vp.push_back(moddefparam("d0","0.7","FULL thickness of shell.","Rsun")); vp.push_back(moddefparam("rb","2.55","Dist to bottom of structure.","Rsun")); vp.push_back(moddefparam("alpha","30.","Angle between axis and foot.","deg")); vp.push_back(moddefparam("rf","10","Dist junction line-circle.","Rsun")); vp.push_back(moddefparam("ratio","0.2","ratio of tube radius to height","Rsun")); return; } // -- density 34 // Test of a simple density rope float CModel34::Density(const Cvec &v) { //float nehomo=density20(v,pparam,temperature); //float returnne=nehomo; float returnne=0; float x=v.v[0],y=v.v[1],z=v.v[2],r=v.norm(); if (r <= 1.01) return 0.; float rbidon,phi,theta; cvcoord(x,y,z,&rbidon,&phi,&theta); // -- Oyz is plane of symmetry x=fabs(x); // -- deal only with first quadran if (y < 0) return returnne; float alpha=pparam[0]; //30.*DTOR; // angle between axis and foot [rad] float rf=pparam[1]; //3.; // dist junction line-circle [Rsun] float tubethick=pparam[2]; //.5; // thickness radius of the rope [Rsun] float constneout=pparam[3]; //1e7; // output constant density float skinthick=pparam[4]; //0.1 // thickness of the skin // -- center of the circle depending on the alpha and rf param // rs : dist solar center - shell skeleton circle center float rs=rf*(cos(alpha)+tan(alpha)*sin(alpha)); // rc : radius shell skeleton circle float rc=rf*tan(alpha); // -- on which side are we ? foot or shell ? float testside=y+x*tan(alpha)-rs; if (testside < 0) { // -- foot side // dist point - foot float r=psldist(Cvec(0,0,0),Cvec(sin(alpha),cos(alpha),0),Cvec(x,y,z)); // constant density for the moment if (r < tubethick && r >= skinthick ) return constneout+returnne; else return returnne; } else { // -- shell side // angular position float rrr,beta,ttt; cvcoord(x,y-rs,0,&rrr,&beta,&ttt); // dist point - circle Cvec vpc(x-rc*cos(beta),y-(rs+rc*sin(beta)),z); // constant density for the moment if (vpc.norm() < tubethick && vpc.norm() >= skinthick) return constneout+returnne; else return returnne; } return returnne; } void CModel34::initParam(float *pparam) { this->pparam=pparam; } void CModel34::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0; vp.push_back(moddefparam("","Simple density rope.","","")); vp.push_back(moddefparam("alpha","0.52","angle between axis and foot","rad")); vp.push_back(moddefparam("rf","3","dist junction line-circle","Rsun")); vp.push_back(moddefparam("tubethick","0.5","thickness radius of the rope","Rsun")); vp.push_back(moddefparam("constneout","1e7","Constant density","cm^-3")); vp.push_back(moddefparam("skinthick","0.1","thickness of the skin","Rsun")); return; } // -- density 35 // Density and temperature cubes from S.Gibson and B.Low CME progam model. // float CModel35::Density(const Cvec &v,float &temperature) { float r,phi,theta; cvcoord(v.v[0],v.v[1],v.v[2],&r,&phi,&theta); float NeighborDensCube[2][2][2]; float NeighborTempCube[2][2][2]; // ---- get the model parameters and point cubes first element float *rco,*phico,*thetaco,*denscube,*tempcube; rco=&pparam[0]; phico=&pparam[3]; thetaco=&pparam[6]; denscube=&pparam[9]; tempcube=&pparam[9+((int) rco[2])*((int) phico[2])*((int) thetaco[2])]; // ---- find the nearest neighbors for tri linear interpolation int mr=(int)(-rco[0]/ rco[1] + r / rco[1] + 1); if (mr >= rco[2] || mr <= 0) return 0.; int mtheta12[2]; int mtheta=(int)(- thetaco[0]/ thetaco[1] + theta / thetaco[1] + 1); if (mtheta >= thetaco[2]) { mtheta12[0]=(int) thetaco[2]-1; mtheta12[1]=mtheta12[0]; } else if (mtheta <= 0) { mtheta12[0]=0; mtheta12[1]=0; } else { mtheta12[0]=mtheta-1; mtheta12[1]=mtheta; } int mphi12[2]; int mphi=(int)(- phico[0]/ phico[1] + phi / phico[1] + 1); if (mphi >= phico[2]) { mphi12[0]=(int) phico[2]-1; mphi12[1]=0; } else if (mphi <= 0) { mphi12[0]=0; mphi12[1]=1; } else { mphi12[0]=mphi-1; mphi12[1]=mphi; } // ---- fill in the nearest neighbor cube for interpolation mr--; mtheta--; mphi--; int xxx,yyy,zzz; for (int ii=0;ii<2;ii++) { xxx=(mr+ii); for (int jj=0;jj<2;jj++) { yyy=mphi12[jj]*((int)(rco[2])); for (int kk=0;kk<2;kk++) { zzz=mtheta12[kk]*((int)(rco[2]) * ((int)(phico[2]))); NeighborDensCube[ii][jj][kk]=*(denscube+(xxx+yyy+zzz)); NeighborTempCube[ii][jj][kk]=*(tempcube+(xxx+yyy+zzz)); } } } // ---- compute intervals for interpolation float rint[2]={rco[0] + rco[1] * mr , rco[0] + rco[1] * (mr+1)}; float phiint[2]={phico[0] + phico[1] * mphi , phico[0] + phico[1] * (mphi+1)}; float thetaint[2]={thetaco[0] + thetaco[1]*mtheta, thetaco[0] + thetaco[1]*(mtheta+1)}; // ---- perform trilinear interpolation float neinterp=trilininterp(r,phi,theta,rint,phiint,thetaint,NeighborDensCube); temperature=trilininterp(r,phi,theta,rint,phiint,thetaint,NeighborTempCube); return neinterp; } void CModel35::initParam(float *pparam) { this->pparam=pparam; } void CModel35::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0; vp.push_back(moddefparam("","Density and temperature cubes from S.Gibson and B.Low CME progam model.","","")); vp.push_back(moddefparam("filedens","''","Density file from gibsoncmewrapper.pro","")); vp.push_back(moddefparam("filetemp","''","Temperature file from gibsoncmewrapper.pro","")); vp.push_back(moddefparam("modparam","","","")); vp.push_back(moddefparam("","load_gibsondens,$filedens,$filetemp,modparam","Load the density and temperature from S.Gibson's model.","")); return; } // -- density 36 // Streamer belt simulation with source surface field map: float CModel36::Density(const Cvec &v) { float x=v.v[0],y=v.v[1],z=v.v[2],r=v.norm(); if (r <= 1.05) return 0.; // -- now deal with orthoradial shape float rbidon,phi,theta; cvcoord(x,y,z,&rbidon,&phi,&theta); // -- param for the neutral sheet map int sang=(int) pparam[0],slat=(int) pparam[1]; // ---- get the distance from the point // to the nearest neighbors // -- nearest neighbor seeking range int srang=30,srlat=30; // in pix // -- point to the neutral sheet map float *pnsheetmap; pnsheetmap=pparam+2; float thetannpos,phinnpos,val;//dist; float thetamg; float nel; int nnok=wherennsmoothed(pnsheetmap,sang,slat,srang,srlat,(phi*RADEG),(theta*RADEG),&phinnpos,&thetannpos,&thetamg,&val); if (nnok == 0) thetamg=DTOR*40; else thetamg*=DTOR; // if (nnok == 0) thetamg=DTOR*((float) srang * 360 / (float) sang); // else thetamg*=DTOR; 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; 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]); } float expo=pow(thetamg/(w*DTOR),2); if (expo > 1E1) nel=np; else nel=np+(ncs-np)*exp(-expo); return nel; } void CModel36::initParam(float *pparam) { this->pparam=pparam; } void CModel36::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=2; vp.push_back(moddefparam("","Streamer belt simulation with source surface field map, user resolution, Guhata. orthoradial.","","")); 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; } // -- density 37 // Streamer belt simulation with source surface field map: float CModel37::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[1],slat=(int) pparam[2]; // -- point to the neutral sheet map float *pnsheetmap; pnsheetmap=pparam+3; float dist,val; // ---- get the distance from the point // here it's given by the magnetic field int nnok=getPosOnSSMap(pnsheetmap,sang,slat,(phi*RADEG),(theta*RADEG),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; } float dcross=w0/(2*u0); // ---- multiply the mag field by a coeff theta=dist*DTOR*pparam[0]; // ca c'est tres bourrin ! 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 CModel37::initParam(float *pparam) { this->pparam=pparam; } void CModel37::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=2; vp.push_back(moddefparam("","Streamer belt simulation with source surface field map, distance from neutral line is given by the magnetic field value.","","")); vp.push_back(moddefparam("sthick","1.","streamer belt thickness coefficient","")); 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($sthick,360,181,nsheetmap,n_elements(nsheetmap))","Format the array of parameters.","")); return; } // -- density 38 // Slab model // see 11 May 2004 // see 14 June 2004 // see 31 Jan 2005 // float CModel38::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); } // ---- FO modulation profile // the *pmodul should contain 120 samples // representing the angles from 60 to 120 deg., // with a range dynamic from 0 to 1 // float maxprofpos=pparam[5]; float aaa=(maxprofpos*(sechang-PI/2))/(2*sechang); float bbb=maxprofpos/(2*sechang); int pos=(int) (aaa+bbb*phi); if (pos < 0) pos=0; if (pos > maxprofpos) pos=(int) maxprofpos; nel*=*(pparam+pos+6); // -- the profile starts at the pos 6 of pparam return nel; } void CModel38::initParam(float *pparam) { this->pparam=pparam; } void CModel38::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { // - pparam[0..3] : coeff of the Ne polynomial // - pparam[4] : slab half angle // - pparam[5] : subscript max of the modulation profile (size - 1) // - pparam[6] : first element of the modulation profile flagcase=0; vp.push_back(moddefparam("","Streamer slab model with parameters useful for fitting.","","")); vp.push_back(moddefparam("c","[1.34e4,1.15e6,-6.022e6,5.577e7]","Coeff of the Ne polynomial.","")); vp.push_back(moddefparam("alpha","0.52","slab half angle","rad")); vp.push_back(moddefparam("smodprof","119","subscript max of the modulation profile (size - 1)","")); vp.push_back(moddefparam("pm","","","")); vp.push_back(moddefparam("","pm=[$c,$alpha,$smodprof,(cos(findgen($smodprof +1)/2)+1.1)/2.1]","","")); return; } // -- density 39 // Test of a simple helix float CModel39::Density(const Cvec &v) { // pos projection point on x axis Cvec op=orthoproj(Cvec(0,0,0),Cvec(1,0,0),v); // distance projection - sun center float dpc=op.norm(); // point of the helix float aa=2,bb=1.; Cvec ph=Cvec(dpc*bb,aa*cos(dpc*0.8),aa*sin(dpc*0.8)); // dist helix center - point Cvec vr=ph-v; if (vr.norm() < 0.5) return 1e7; else return 0; //float nehomo=density20(v,pparam,temperature); //float returnne=nehomo; float returnne=0; float x=v.v[0],y=v.v[1],z=v.v[2],r=v.norm(); if (r <= 1.01) return 0.; float rbidon,phi,theta; cvcoord(x,y,z,&rbidon,&phi,&theta); // -- Oyz is plane of symmetry x=fabs(x); // -- deal only with first quadran if (y < 0) return returnne; float alpha=pparam[0]; //30.*DTOR; // angle between axis and foot [rad] float rf=pparam[1]; //3.; // dist junction line-circle [Rsun] float tubethick=pparam[2]; //.5; // thickness radius of the rope [Rsun] float constneout=pparam[3]; //1e7; // output constant density // -- center of the circle depending on the alpha and rf param // rs : dist solar center - shell skeleton circle center float rs=rf*(cos(alpha)+tan(alpha)*sin(alpha)); // rc : radius shell skeleton circle float rc=rf*tan(alpha); // -- on which side are we ? foot or shell ? float testside=y+x*tan(alpha)-rs; if (testside < 0) { // -- foot side // point on the foot Cvec op=orthoproj(Cvec(0,0,0),Cvec(sin(alpha),cos(alpha),0),Cvec(x,y,z)); // helix position float aa=0.2,omega=0.1; // param of the helix Cvec hp=Cvec(aa*cos(op.norm()*omega),aa*sin(op.norm()*omega),0)+op; // dist point - helix Cvec vr=Cvec(x,y,z)-hp; // constant density for the moment if (vr.norm() < tubethick) return constneout+returnne; else return returnne; } else { // -- shell side // angular position float rrr,beta,ttt; cvcoord(x,y-rs,0,&rrr,&beta,&ttt); // dist point - circle Cvec vpc(x-rc*cos(beta),y-(rs+rc*sin(beta)),z); // constant density for the moment if (vpc.norm() < tubethick) return constneout+returnne; else return returnne; } return returnne; } void CModel39::initParam(float *pparam) { this->pparam=pparam; } void CModel39::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0; vp.push_back(moddefparam("","Test of a simple helix.","","")); vp.push_back(moddefparam("alpha","0.52","angle between axis and foot","rad")); vp.push_back(moddefparam("rf","3.","dist junction line-circle","Rsun")); vp.push_back(moddefparam("tubethick","0.5","thickness radius of the rope","Rsun")); vp.push_back(moddefparam("constneout","1e7","output constant density","cm^-3")); return; } // -- density 40 // Streamer belt simulation with source surface field map // float CModel40::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[0],slat=(int) pparam[1]; // -- point to the neutral sheet map float *pnsheetmap; pnsheetmap=pparam+3; float dist,val; // ---- get the distance from the point // here it's given by the magnetic field int nnok=getPosOnSSMap(pnsheetmap,sang,slat,(phi*RADEG),(theta*RADEG),dist,val); if (nnok == 0) return 0.; // ---- get the height of the passed Ne layer // float layerheigth=pparam[2]; // ---- compute Saito model (equatorial) // float c1=1.36e6; // float d1=2.14; // float c2=1.68e8; // float d2=6.13; // ---- compute ratio for scaling //if (val < 7e5 ) return 0.; //float Neratio=val/(c1*pow(layerheigth,-d1)+c2*pow(layerheigth,-d2)); // ---- compute Ne and scale //float nel=Neratio*(c1*pow(r,-d1)+c2*pow(r,-d2)); //float nel=(c1*pow(r,-d1)+c2*pow(r,-d2)); 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; } 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); } // ---- modulate by the layer value return val*nel; } void CModel40::initParam(float *pparam) { this->pparam=pparam; } void CModel40::dumpDefaultParamForIDL(std::vector& vp,int& flagcase) { flagcase=0; //! - pparam[0] : longitude size in pix //! - pparam[1] : latitude size in pix //! - pparam[2] : height of the layer in Rsun //! - pparam[3] : Ne layer map (lon,lat) vp.push_back(moddefparam("","Streamer belt simulation with source surface field map, Saito radial scaled on the Ne layer map passed in parameter.","","")); vp.push_back(moddefparam("pfssmapfile","'WSOpsffR250CR2012.fts'","!!! ACTUALLY THE MAP SHOULD BE NE NOT A MAG FIELD !!! THIS IS NOT A VALID EXAMPLE !!!","")); vp.push_back(moddefparam("heightpfss","2.5","","")); vp.push_back(moddefparam("mp","","","")); vp.push_back(moddefparam("","if not file_exist($pfssmapfile) then begin","","")); vp.push_back(moddefparam("","rtinitenv","","")); vp.push_back(moddefparam("","progpath=getenv('RT_PATH')+get_delim()","","")); vp.push_back(moddefparam("","endif else progpath=''","","")); vp.push_back(moddefparam("","map=readfits(progpath + $pfssmapfile)","Load the pfss map","")); vp.push_back(moddefparam("","mp=[(size(map,/dim))[0],(size(map,/dim))[1],$heightpfss,reform(map,n_elements(map))]","","")); return; } /* * $Log: models31to40.cpp,v $ * Revision 1.3 2009/02/09 20:51:15 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 * */