// $Id: voxinteg.h,v 1.2 2007/05/14 17:19:41 thernis Exp $ { // ---- compute the vertices pos on the image voxelvert2iijj(vij,voxel.vertobs,xstep,ystep,icntr,jcntr); // ---- find the extrema of pix pos in the image vector::const_iterator it= min_element(vij[0].begin(),vij[0].end()); float iminf=*it; int imin=(int) (*it +0.5); // smallest element in I it = min_element(vij[1].begin(),vij[1].end()); float jminf=*it; int jmin=(int) (*it +0.5); // smallest element in J it = max_element(vij[0].begin(),vij[0].end()); float imaxf=*it; int imax=(int) (*it +1); // largest element in I it = max_element(vij[1].begin(),vij[1].end()); float jmaxf=*it; int jmax=(int) (*it +1); // largest element in J // -- check if not out of image bounds if (imin>=*is || jmin>=*js || imax<0 || jmax<0) continue; if (imin < 0) { imin=0; iminf=-0.5; } if (imax >= *is) { imax=*is-1; imaxf=float(imax)+0.5; } if (jmin < 0) { jmin=0; jminf=-0.5; } if (jmax >= *js) { jmax=*js-1; jmaxf=float(jmax)+0.5; } /* // ---- Test if voxel under resolution of 1 pix // -0.5 for shift to the edge of the pix: 0,0 is the cntr of the pix if ( imaxf-floor(iminf+0.5)+0.5<1 && jmaxf-floor(jminf+0.5)+0.5<1) { countunder++; // ---- just use the vox center pos and total Ne int ip=int(iminf+0.5); int jp=int(jminf+0.5); if (ip < 0 || ip >= *is || jp < 0 || jp >= *js) continue; int impos=jp * *is + ip; // -- compute Ne netot[impos]+=voxel.dens*voxel.volume; // -- loop if only Ne is requested if (*flagneonly) continue; // -- impact parameter of the voxel center Cvec Vvoxcenterabs=obs.ui*voxel.voxcenterobs; float rho=psldist(obs.o,Vvoxcenterabs,abs.o); // -- skip integration if LOS within the occulter if (rho < *poccrad) continue; // -- dist vox cntr - sun cntr float r=(Vvoxcenterabs+obs.o).mag(); // -- integration only in front of the Sun if (!(rho > LIMBLIMIT || (r > LIMBLIMIT && Vvoxcenterabs[2] < 0))) continue; float lengththru=0; voxel.GetThruLengthCenter(lengththru,0); // -- Compute Thomson scattering geometric factors and then the brightness float totterm,polterm; getThomsonGeomFactor(r, rho, totterm,polterm); btot[impos]+=constfactor*lengththru*totterm*voxel.dens; bpol[impos]+=constfactor*lengththru*polterm*voxel.dens; continue; } */ /* // -- Test if voxel under resolution. // If it falls into that case it's that the pix is in // between two pix. // Compute barycenter of vox projection and spread the density in a // maxrespix boxcar. // -0.5 for shift to the edge of the pix: 0,0 is the cntr of the pix // for the moment use a 2x2 boxcar if ( false && imaxf-iminf<=1 && jmaxf-jminf<=1) { countunder2++; // ---- just use the vox center pos and total Ne // -- compute projection of the vox center on the image float vij[2]; voxelvert2iijj(vij,voxel.voxcenterobs,xstep,ystep,icntr,jcntr); int ip=int(vij[0]); int jp=int(vij[1]); int impos=jp * *is + ip; if (impos < maxnbpix) { float xw=(vij[0]-float(ip)); float yw=(vij[1]-float(jp)); // -- compute Ne netot[impos]+=voxel.dens*((1-xw)*(1-yw))*voxel.volume; if ((imin+1) < *is) netot[impos+1]+=voxel.dens*(xw*(1-yw))*voxel.volume; if (imin+1 < *is && jmin+1 < *js) netot[impos+1+ *is]+=voxel.dens*(xw*yw)*voxel.volume; if (jmin+1 < *js) netot[impos+ *is]+=voxel.dens*((1-xw)*yw)*voxel.volume; } // -- loop if only Ne is requested if (*flagneonly) continue; // -- impact parameter of the voxel center Cvec Vvoxcenterabs=obs.ui*voxel.voxcenterobs; float rho=psldist(obs.o,Vvoxcenterabs,abs.o); // -- skip integration if LOS within the occulter if (rho < *poccrad) continue; // -- dist vox cntr - sun cntr float r=(Vvoxcenterabs+obs.o).mag(); // -- integration only in front of the Sun if (!(rho > LIMBLIMIT || (r > LIMBLIMIT && Vvoxcenterabs[2] < 0))) continue; float lengththru=0; voxel.GetThruLengthCenter(lengththru,0); // -- Compute Thomson scattering geometric factors and then the brightness float totterm,polterm; getThomsonGeomFactor(r, rho, totterm,polterm); if (impos < maxnbpix) { float xw=(vij[0]-float(ip)); float yw=(vij[1]-float(jp)); // -- compute btot and bpol btot[impos]+=constfactor*lengththru*totterm*voxel.dens*((1-xw)*(1-yw)); bpol[impos]+=constfactor*lengththru*polterm*voxel.dens*((1-xw)*(1-yw)); rrr[impos]=lengththru; if ((imin+1) < *is) { btot[impos+1]+=constfactor*lengththru*totterm*voxel.dens*(xw*(1-yw)); bpol[impos+1]+=constfactor*lengththru*polterm*voxel.dens*(xw*(1-yw)); } if (imin+1 < *is && jmin+1 < *js) { btot[impos+1+ *is]+=constfactor*lengththru*totterm*voxel.dens*(xw*yw); bpol[impos+1+ *is]+=constfactor*lengththru*polterm*voxel.dens*(xw*yw); } if (jmin+1 < *js) { btot[impos+ *is]+=constfactor*lengththru*totterm*voxel.dens*((1-xw)*yw); bpol[impos+ *is]+=constfactor*lengththru*polterm*voxel.dens*((1-xw)*yw); } } continue; } */ // -- dist obs to voxcenter float distobs2vox=voxel.voxcenterobs.mag(); // -- surface area corresponding to the pixel square cone at the distance of the observer float voxpixedge=2*distobs2vox*tan(*fovpix/2); //*RSUN_CM float voxpixsurf=voxpixedge*voxpixedge; // ---- Scan the LOSes going through the voxel Cvec Vcntrvoxel; float lengththru; Cvec Vsuncntrvoxel=voxel.base.u*(voxel.base.o-abs.o); // pos of the sun cntr in the voxel coord sys countbig++; float totterm,polterm; // thomson scattering geometrical factors for (int jj=jmin;jj<=jmax;jj++) { int jpos=jj * *is; for (int ii=imin;ii<=imax;ii++) { int impos=ii + jpos; if (vlosimobs[impos].IsNull()) { vlosimobs[impos]=iijj2losobs((float) ii,(float) jj,xstart,ystart,xstep,ystep); } int outvoxelflag=voxel.GetThruLength(vlosimobs[impos],Vcntrvoxel,lengththru,0); // -- loop if out LOS out of the voxel if (outvoxelflag) continue; // -- compute Ne netot[impos]+=voxel.dens*lengththru*voxpixsurf; // -- loop if only Ne is requested if (*flagneonly) continue; if (vlosimabs[impos].IsNull()) { vlosimabs[impos]=obs.ui * vlosimobs[impos]; rho[impos]=psldist(obs.o,vlosimabs[impos],abs.o); } // -- skip integration if LOS within the occulter if (rho[impos] < *poccrad) continue; // -- dist vox cntr - sun cntr float r=(Vsuncntrvoxel+Vcntrvoxel).mag(); // -- integration only in front of the Sun if (!(rho[impos] > LIMBLIMIT || (r > LIMBLIMIT && vlosimabs[impos][2] < 0))) continue; // -- Compute Thomson scattering geometric factors and then the brightness getThomsonGeomFactor(r, rho[impos], totterm,polterm); btot[impos]+=constfactor*lengththru*totterm*voxel.dens; bpol[impos]+=constfactor*lengththru*polterm*voxel.dens; } } } /* * $Log: voxinteg.h,v $ * Revision 1.2 2007/05/14 17:19:41 thernis * Add CVS id and log in all files * */