/*************************************************************************** * raytrace.h * * $Id: raytrace.h,v 1.16 2010/09/01 19:51:18 thernis Exp $ ****************************************************************************/ /** \file raytrace.h * \brief Contains all the wrappers for IDL */ #ifndef RAYTRACE_H #define RAYTRACE_H #ifdef HAVE_CONFIG_H #include #endif #include "constant.h" #include "cuvemission.h" #include "Cvec.h" #include "Clos.h" #include "CModelBase.h" #include "Cdvoxel.h" #include "rtmiscfunc.h" #include "ModelPosition.h" //! Raytracing program parameters, for parameter passing to raytracemain struct rtparam { int *pmodelid; //!< Model number int *pis; //!< X size of the image int *pjs; //!< Y size of the image float *pfovpix; //!< Pixel resolution [rad] float *pobspos; //!< Observer position [Rsun] float *pobsang; //!< Observer camera orientation [rad] float *pnepos; //!< Electron density center positon [Rsun] float *pneang; //!< Electron density orientation [rad] int *plosnbp; //!< LOS number of integration steps float *plosrange; //!< LOS integration range [Rsun] float *pmparam; //!< Model parameter array pointer int *pquiet; //!< Flag: Disable display of program progression float *pbtot; //!< Total brightness image pointer float *pbpol; //!< Polarized brightness image pointer float *pnetot; //!< Integrated electron density image pointer float *pcrpix; //!< Position of the Sun center on the image float *prhoim; //!< Image giving the impact parameter float *pmmlon; //!< Min and max longitude float *pmmlat; //!< Min and max latitude float *prrr; //!< Projected radius on the Earth plane of the sky int *pneonly; //!< Flag to set for Ne only calculation int *proi; //!< Region of interest map int *ppofinteg; //!< Force integration centered on a plane of integration float *ppoiang; //!< Plane of integration angle orientation float *Hlonlat; //!< Heliographic lon and lat of disk center in rad float *poccrad; //!< Radius of the occulter. float *padapthres; //!< Adaptative Simpson's integration difference threshold. No Simpson integration if = 0. int *pmaxsubdiv; //!< Adaptative Simpson's integration maximun number of interval subdivision recursion. float *plimbdark; //!< User limb darkening: default 0.58 float *protmat; //!< output final Ne rotation matrix float *obslonlat; //!< observer lon lat and height in carrington coordinate int *obslonlatflag; //!< set to 1 if user defined obslonlat: then obspos is ignored int *projtypecode; //!< type of projection: 1:ARC, 2:TAN, 3:SIN, 4:AZP float *pv2_1; //!< mu parameter for the AZP projection int *pfrontinteg; //!< Integration centered on the observer unsigned int uvinteg; //!< Use UV emission instead of Thomson scattering float disttofracmax; //!< Set to the fraction of max B to evaluate the integration distance to that threshold: the distance is returned in bpol float *pnerotcntr; //!< Electron density center positon within the nps coord system defined by nepos and neang [Rsun] float *pnerotang; //!< Electron density orientation within the nps coord system defined by nepos and neang [rad] float *pnetranslation; //!< Shift of the electron density, in the density coordinate system int *pnerotaxis; //!< Electron density axis ids corresponding to the nerotang rotation angles }; //! Main raytracing routine in cartesian coordinates int raytracemain(rtparam fp); //! Wrapper of raytracemain called by IDL extern "C" int raytracewl(int argc, void **argv); //! Raytracing using WCS coordinate system extern "C" int rtraytracewcs(int argc, void **argv); //! Raytracing program that integrates only a cirular profile extern "C" int rtwlcirc(int argc, void **argv); //! raytracing program that integrates only along a straight line segment profile extern "C" int rtwlseg(int argc, void **argv); //! raytracing for only for a density cube extern "C" int rtdenscube(int argc, void **argv); //! return the carrington position on the solar limb for a given pixel: useful for EUVI extern "C" int rtGetCarPosOnDisk(int argc, void **argv); //! Compute ii,jj pix pos for the cube vertices std::vector > voxelvert2iijj(std::vector &vert,const float &xstep,const float &ystep,const float &icntr,const float &jcntr); //! Compute ii,jj pix pos for the cube vertices std::vector voxelvert2iijj(Cvec &vert,const float &xstep,const float &ystep,const float &icntr,const float &jcntr); //! Convert LOS orientation in obs to pix coordinates void losobs2xxyy(const Cvec &Vlosobs,float &xx,float &yy); void losobs2xxyydebug(const Cvec &Vlosobs,float &xx,float &yy); //! Convert pix coordinates to LOS orientation in obs inline void xxyy2losobs(const float &xx,const float &yy,float &rrr,float &alpha) { rrr=sqrt(xx*xx+yy*yy); alpha=0; if (yy != 0) { alpha=atan(xx/yy); if (yy < 0) alpha=alpha+PI; } } //! Convert pix coordinates to LOS orientation in obs inline Cvec xxyy2losobs(float const &xx,float const &yy) { float rrr,alpha; xxyy2losobs(xx,yy,rrr,alpha); return Cvec(sin(rrr)*sin(alpha),sin(rrr)*cos(alpha),cos(rrr)); } //! Convert pix coordinates to LOS orientation in obs inline void xxyy2losobs(const float &i,const float &j,float &rrr,float &alpha,const float *crpix,const float *fovpix) { float x=i-crpix[0]; float y=j-crpix[1]; rrr=*fovpix*sqrt(x*x+y*y); alpha=atan2(y,x); } //! Convert pix coordinates to LOS orientation in obs inline void xxyy2losobs(const float &i,const float &j,float &rrr,float &alpha,const float *crpix,const float *fovpix,const float *pc) { // ---- see Calabretta & Greisen, A&A, 2002 float ic=i-crpix[0]; float jc=j-crpix[1]; float x=pc[0]*ic+pc[2]*jc; float y=pc[1]*ic+pc[3]*jc; rrr=*fovpix*sqrt(x*x+y*y); alpha=atan2(y,x); } //! Convert LOS (in obs coord) into pix position inline void losobs2xxyy(const Cvec &Vlosobs,const float *crpix,const float *fovpix,const float *pci,const int &projtypecode,const float &pv2_1,float &i, float &j) { float phi=atan2( Vlosobs.v[0], Vlosobs.v[1]); // float theta0=acos( Vlosobs.v[2] / Vlosobs.norm()); float theta0=atan2(sqrt(Vlosobs.v[0]*Vlosobs.v[0]+Vlosobs.v[1]*Vlosobs.v[1]),Vlosobs.v[2]); float theta=applyprojection(theta0,projtypecode,pv2_1); // see 10 Jul 2007 float r=theta / (*fovpix); float x=r*cos(phi); float y=r*sin(phi); float ic=pci[0]*x+pci[2]*y; float jc=pci[1]*x+pci[3]*y; i=ic+crpix[0]; j=jc+crpix[1]; } //! Convert pix coordinates to LOS orientation in obs inline Cvec xxyy2losobs(const float &i,const float &j,float &rrr,float &alpha,const float *crpix,const float *fovpix,const float *pc,const int &projtypecode,const float &pv2_1) { // ---- see Calabretta & Greisen, A&A, 2002 float ic=i-crpix[0]; float jc=j-crpix[1]; float x=pc[0]*ic+pc[2]*jc; float y=pc[1]*ic+pc[3]*jc; rrr=*fovpix*sqrt(x*x+y*y); alpha=atan2(y,x); rrr=applyinverseprojection(rrr,projtypecode,pv2_1); return Cvec(sin(rrr)*sin(alpha),sin(rrr)*cos(alpha),cos(rrr)); } //! returns the determinant and the two roots (if det >= 0) of a 2nd deg polynomial. a0+a1x+a2x^2 inline float get2ndDegPolyRoots(const float &a0,const float &a1,const float &a2,float &root1,float &root2) { float det=a1*a1-4.*a0*a2; if (det < 0.) return det; float sqrtdet=sqrt(det); root1=-(a1+sqrtdet)/(2.*a2); root2=(-a1+sqrtdet)/(2.*a2); return det; } //! Compute impact parameter and LOS direction inline void getimpactandlos(const float &xx,const float &yy,const Cbasis &obs,const Cbasis &abs,Cvec &vlosabs,float &rho) { float rrr; float alpha=0; xxyy2losobs(xx,yy,rrr,alpha); vlosabs=obs.ui * Cvec(sin(rrr)*sin(alpha),sin(rrr)*cos(alpha),cos(rrr)); // 3 // -- compute rho : dist LOS - Sun cntr: impact parameter rho=psldist(obs.o,vlosabs,abs.o); } //! Compute impact parameter and LOS direction inline void getimpactandlos(const float &xx,const float &yy,const Cbasis &obs,const Cbasis &abs,const int projtypecode,const float pv2_1,Cvec &vlosabs,float &rho) { float rrr; float alpha=0; // ---- compute LOS according to pixel position xxyy2losobs(xx,yy,rrr,alpha); // -- apply inverse WCS projection rrr=applyinverseprojection(rrr,projtypecode,pv2_1); vlosabs=obs.ui * Cvec(sin(rrr)*sin(alpha),sin(rrr)*cos(alpha),cos(rrr)); // 3 // -- compute rho : dist LOS - Sun cntr: impact parameter rho=psldist(obs.o,vlosabs,abs.o); } //! Compute impact parameter and LOS direction inline void getimpactandlos(const float &i,const float &j,const Cbasis &obs,const Cbasis &abs,const int projtypecode,const float pv2_1,Cvec &vlosabs,float &rho,const float *crpix,const float *fovpix,const float *pc) { float rrr; float alpha=0; // ---- get the los according to the pixel position xxyy2losobs(i,j,rrr,alpha,crpix,fovpix,pc); // -- apply inverse WCS projection rrr=applyinverseprojection(rrr,projtypecode,pv2_1); vlosabs=obs.ui * Cvec(sin(rrr)*sin(alpha),sin(rrr)*cos(alpha),cos(rrr)); // -- compute rho : dist LOS - Sun cntr: impact parameter rho=psldist(obs.o,vlosabs,abs.o); } //! Compute impact parameter and LOS direction inline void getimpactandlos(const float &i,const float &j,const Cbasis &obs,const Cbasis &abs,const int projtypecode,const float pv2_1,Cvec &vlosabs,float &rho,float &rrr,float &alpha,Cvec &vlosobs,const float *crpix,const float *fovpix,const float *pc) { // ---- get the los according to the pixel position xxyy2losobs(i,j,rrr,alpha,crpix,fovpix,pc); printvar(*fovpix); printvar(pc[0]); printvar(rrr); printvar(alpha); // -- apply inverse WCS projection rrr=applyinverseprojection(rrr,projtypecode,pv2_1); vlosobs=Cvec(sin(rrr)*sin(alpha),sin(rrr)*cos(alpha),cos(rrr)); vlosabs=obs.ui * vlosobs; // -- compute rho : dist LOS - Sun cntr: impact parameter rho=psldist(obs.o,vlosabs,abs.o); printvar(crpix[0]); printvar(crpix[1]); printvar(i); printvar(j); } //! Compute the LOS corresponding to the image pixel position inline Cvec iijj2losobs(const float &i,const float &j,const float &xstart,const float &ystart,const float &xstep,const float &ystep) { float xx=xstart+j*xstep; float yy=ystart+i*ystep; return xxyy2losobs(xx,yy); } //! Compute ii,jj pix pos for the voxel vertices inline void voxelvert2iijj(std::vector > &vij,std::vector &vert,const float &xstep,const float &ystep,const float &icntr,const float &jcntr) { float xx,yy; for(unsigned int i=0;i sqrt(obs.o.mag()*obs.o.mag()-rho*rho) && pscal(vrho-obs.o,vs-obs.o) >= 0)) continue; // if (rho > LIMBLIMIT || (r > LIMBLIMIT && vs.v[2] > 0.)) /******************************************/ /* MODEL CALL HERE */ /******************************************/ // -- compute Ne //float ner=pmod->Density((nps.u * (vs - nps.o))-nps.translation); float ner=pmod->Density(ChangetoDensityCoord(modpos,vs)); if (ner <= 1e-1) continue; // total density *posne +=ner; // -- skip brightness calculation if requested if (flagneonly > 0) continue; /******************************************/ // -- compute geometrical functions float sinomega=RSUN/r; float sinsquareomega=sinomega*sinomega; float cossquareomega=1-sinsquareomega; float cosomega=sqrt(cossquareomega); float logterm=log((1.+sinomega)/cosomega); float a=cosomega*sinsquareomega; float b=-1./8.*(1.-3.*sinsquareomega-cossquareomega*((1.+3.*sinsquareomega)/sinomega)*logterm); float c=(4./3.)-cosomega-(cosomega*cossquareomega)/3.; float d=(1./8.)*(5.+sinsquareomega-cossquareomega*((5.-sinsquareomega)/sinomega)*logterm); // ---- sum in the pixel float rhooverr=rho/r; // the polarized brightness float polterm=(a+u*(b-a))*rhooverr*rhooverr; *posbpol +=ner*polterm; // the total brightness *posbtot +=ner*(2*(c+u*(d-c))-polterm); } // multiply by the integral constant factor *posbtot *=constfactor*los.ds; *posbpol *=constfactor*los.ds; *posne *=RSUN_CM*los.ds; } //! Check if a chunk is not within or behind the sun inline bool checkchunkwbsun(float rho,float r0, float r2,float vs0z,float vs2z) { return (!((rho > LIMBLIMIT) || ((r0 > LIMBLIMIT && vs0z < 0) && (r2 > LIMBLIMIT && vs2z < 0)))); } //! Calculate and return the Thomson scattering geometric functions inline void getThomsonGeomFactor(const float &r,const float &rho,float &totterm,float &polterm) { float sinomega=RSUN/r; float sinsquareomega=sinomega*sinomega; float cossquareomega=1-sinsquareomega; float cosomega=sqrt(cossquareomega); float logterm=log((1.+sinomega)/cosomega); float a=cosomega*sinsquareomega; float b=-1./8.*(1.-3.*sinsquareomega-cossquareomega*((1.+3.*sinsquareomega)/sinomega)*logterm); float c=(4./3.)-cosomega-(cosomega*cossquareomega)/3.; float d=(1./8.)*(5.+sinsquareomega-cossquareomega*((5.-sinsquareomega)/sinomega)*logterm); float rhooverr=rho/r; // Polarized brightness term polterm=(a+U*(b-a))*rhooverr*rhooverr; // Total brightness term totterm=(2*(c+U*(d-c))-polterm); } //! Integration with adaptative step inline void adaptabacksim(float rvar,Cvec vs,float f0ne,float f1ne,float f2ne,float losstep,Cvec vlosstep,int depthcount,float &sum,float &sbtot,float &sbpol,CModelBase* pmod,float* pmodparam,ModelPosition modpos,int &depth,int maxsubdiv,float rho) { float f01,f12,sum01=0.,sum12=0.,sbtot01=0.,sbtot12=0., sbpol01=0.,sbpol12=0.; float tott[3],polt[3]; // -- give up if too many recusions if (depthcount > maxsubdiv) { sum=losstep*(f0ne+4*f1ne+f2ne)/3; // -- get Thomson scattering geometric factors getThomsonGeomFactor(vs.norm(),rho,tott[0],polt[0]); getThomsonGeomFactor((vs+vlosstep).norm(),rho,tott[1],polt[1]); getThomsonGeomFactor((vs+(vlosstep*2)).norm(),rho,tott[2],polt[2]); sbtot=losstep*(f0ne*tott[0]+4*f1ne*tott[1]+f2ne*tott[2])/3; sbpol=losstep*(f0ne*polt[0]+4*f1ne*polt[1]+f2ne*polt[2])/3; //cout << "---- Max Depth : " << depthcount << endl; depth=depthcount; return; } // -- check if difference too big if (fabs(f2ne-f0ne) > rvar) { // -- divide the interval if too big // -- left // f01=pmod->Density(nps.u * ((vs+(vlosstep*0.5)) + nps.o)); f01=pmod->Density(ChangetoDensityCoord(modpos,vs+(vlosstep*0.5))); adaptabacksim(rvar,vs,f0ne,f01,f1ne,losstep*0.5,vlosstep*0.5,depthcount+1,sum01,sbtot01,sbpol01,pmod,pmodparam,modpos,depth,maxsubdiv,rho); // -- divide the interval if too big // -- right // f12=pmod->Density(nps.u * ((vs+(vlosstep*1.5)) + nps.o)); f12=pmod->Density(ChangetoDensityCoord(modpos,vs+(vlosstep*1.5))); adaptabacksim(rvar,vs+vlosstep,f1ne,f12,f2ne,losstep*0.5,vlosstep*0.5,depthcount+1,sum12,sbtot12,sbpol12,pmod,pmodparam,modpos,depth,maxsubdiv,rho); // -- sum both sides sum=sum01+sum12; sbtot=sbtot01+sbtot12; sbpol=sbpol01+sbpol12; } else { // ---- sum if OK // -- Ne sum=losstep*(f0ne+4*f1ne+f2ne)/3; // -- get Thomson scattering geometric factors getThomsonGeomFactor(vs.norm(),rho,tott[0],polt[0]); getThomsonGeomFactor((vs+vlosstep).norm(),rho,tott[1],polt[1]); getThomsonGeomFactor((vs+(vlosstep*2)).norm(),rho,tott[2],polt[2]); sbtot=losstep*(f0ne*tott[0]+4*f1ne*tott[1]+f2ne*tott[2])/3; sbpol=losstep*(f0ne*polt[0]+4*f1ne*polt[1]+f2ne*polt[2])/3; depth=depthcount; } return; } //! Integration along a line of sight with adaptative step inline void losintegadaptstep(const int &pofinteg, const int &frontinteg, const Cbasis &obs, const Cvec &vlosabs, const Cbasis &abs, Cbasis &poi, Clos &los, const float &rho, CModelBase *pmod, float *pmodparam, const ModelPosition &modpos, const float &u, const float &constfactor, const int &flagneonly, const float &adapthres, const int &maxsubdiv, float *posbtot, float *posbpol, float *posne) { Cvec qlos; calcqlos(pofinteg,frontinteg,obs,vlosabs,abs,poi,qlos); // -- initialize *posbtot=0; *posbpol=0; *posne=0; // -- current position on the LOS float s=los.sstart; Cvec vlosstep=vlosabs * los.ds; Cvec vs=qlos + vlosabs * s; // -- loop bool flagfirstloop=true; float f0ne=0,f1ne,f2ne; for(unsigned int k=0;kDensity(nps.u * (vs + nps.o)); f0ne=pmod->Density(ChangetoDensityCoord(modpos,vs)); flagfirstloop=false; } // -- compute Ne : f1, f2 /* f1ne=pmod->Density(nps.u * ((vs+(vlosstep*0.5)) + nps.o)); f2ne=pmod->Density(nps.u * (vs2 + nps.o));*/ f1ne=pmod->Density(ChangetoDensityCoord(modpos,vs+(vlosstep*0.5))); f2ne=pmod->Density(ChangetoDensityCoord(modpos,vs2)); // -- call backtrack integrator float sumne02=0,stot=0,spol=0; int depth=0; adaptabacksim(adapthres,vs,f0ne,f1ne,f2ne,los.ds*0.5,vlosstep*0.5,1,sumne02,stot,spol,pmod,pmodparam,modpos,depth,maxsubdiv,rho); // -- add density *posne +=sumne02; *posbtot +=stot; *posbpol +=spol; // -- switch last to first point f0ne=f2ne; } // multiply by the integral constant factor *posbtot *=constfactor; *posbpol *=constfactor; *posne *=RSUN_CM; } //! Integration along a line of sight for UV light inline void losintegUV(const int &pofinteg, const int &frontinteg, const Cbasis &obs, const Cvec &vlosabs, const Cbasis &abs, Cbasis &poi, Clos &los, const float &rho, CModelBase *pmod, float *pmodparam, const ModelPosition &modpos, const float &u, const float &constfactor, const int &flagneonly, float *posbtot, float *posbpol, float *posne, CUVEmission *uvemis, const int &uvinteg) { // -- Compute the origin of the LOS, depending what the user specified Cvec qlos; calcqlos(pofinteg,frontinteg,obs,vlosabs,abs,poi,qlos); // -- initialize *posbtot=0; *posbpol=0; *posne=0; // -- current position on the LOS float s=los.sstart; Cvec vlosstep=vlosabs * los.ds; Cvec vs=qlos + vlosabs * s; float temperature; for(unsigned int k=0;k LIMBLIMIT || (r > LIMBLIMIT && vs.v[2] > 0.)) { /******************************************/ /* MODEL CALL HERE */ /******************************************/ // -- compute Ne // float ner=pmod->Density(nps.u * (vs + nps.o),temperature); float ner=pmod->Density(ChangetoDensityCoord(modpos,vs),temperature); float uv=uvemis->calcEmissivity(uvinteg,temperature); // total density *posne +=ner; *posbpol +=uv; *posbtot +=uv*ner*ner; } } // multiply by the integral constant factor *posbtot *=los.ds; *posbpol *=los.ds; *posne *=RSUN_CM*los.ds; } //! Integration along a line of sight inline void losintegdisttofracmax(const int &pofinteg, const int &frontinteg, const Cbasis &obs, const Cvec &vlosabs, const Cbasis &abs, Cbasis &poi, Clos &los, const float &rho, CModelBase *pmod, float *pmodparam, const ModelPosition &modpos, const float &u, const float &constfactor, const int &flagneonly, float *posbtot, float *posbpol, float *posne, float disttofracmax) { // -- Compute the origin of the LOS, depending what the user specified Cvec qlos; calcqlos(pofinteg,frontinteg,obs,vlosabs,abs,poi,qlos); // -- initialize: we return the half distance value in the bpol image *posbpol=0; // -- current position on the LOS float bucket=0.; Cvec vlosstep=vlosabs * los.ds; Cvec vs=qlos + vlosabs * los.sstart; Cvec vrho=orthoproj(obs.o,vlosabs,abs.o); for(unsigned int k=0;k sqrt(obs.o.mag()*obs.o.mag()-rho*rho) && pscal(vrho-obs.o,vs-obs.o) >= 0)) continue; // -- integration only in front of the Sun // if (rho > LIMBLIMIT || (r > LIMBLIMIT && vs.v[2] > 0.)) { /******************************************/ /* MODEL CALL HERE */ /******************************************/ // -- compute Ne float temperature; float ner=pmod->Density(ChangetoDensityCoord(modpos,vs)); if (ner <= 1e-1) continue; /******************************************/ // -- compute geometrical functions float sinomega=RSUN/r; float sinsquareomega=sinomega*sinomega; float cossquareomega=1-sinsquareomega; float cosomega=sqrt(cossquareomega); float logterm=log((1.+sinomega)/cosomega); float a=cosomega*sinsquareomega; float b=-1./8.*(1.-3.*sinsquareomega-cossquareomega*((1.+3.*sinsquareomega)/sinomega)*logterm); float c=(4./3.)-cosomega-(cosomega*cossquareomega)/3.; float d=(1./8.)*(5.+sinsquareomega-cossquareomega*((5.-sinsquareomega)/sinomega)*logterm); // ---- sum in the pixel float rhooverr=rho/r; // the polarized brightness float polterm=(a+u*(b-a))*rhooverr*rhooverr; // the total brightness bucket +=ner*(2*(c+u*(d-c))-polterm)*constfactor*los.ds; if (bucket > (disttofracmax * (*posbtot))) break; *posbpol+=los.ds; // } } } #endif /* * $Log: raytrace.h,v $ * Revision 1.16 2010/09/01 19:51:18 thernis * Implement disttofracmax * * Revision 1.15 2009/04/13 21:03:57 thernis * - Use ModelPosition class. * - Implement extra positioning parameters for the models. * * Revision 1.14 2009/03/06 21:23:46 thernis * Implement neshift * * Revision 1.13 2009/02/23 16:21:15 thernis * Fix calculation of brightness in the UV renderer * * Revision 1.12 2009/02/09 20:51:20 thernis * - Clean up the code * - Change CModel::Density prototype * - Update documentation * - Implement multi-threading using boost thread * - Add new models * * Revision 1.11 2008/09/30 18:33:29 thernis * Replace include of losinteg.h and losintegadaptstep.h by inline functions. * * Revision 1.10 2008/09/23 14:08:13 thernis * - add new testing models, - implement integration in front of the instrument * * Revision 1.9 2007/07/19 19:51:03 thernis * Implement rtGetCarPosOnDisk * * Revision 1.8 2007/07/10 21:15:57 thernis * Implement raytracing of a cloud of points. * * Revision 1.7 2007/05/14 17:19:41 thernis * Add CVS id and log in all files * */