pro for_pinpoint,pB,pol,r3D,tau3D,Density,rpos,distobs,losuse,xx_plus,xx_minus,xx_ground,tau_plus,tau_minus,tau_ground,ulimb,$ dimrad=dimrad,thresdens=thresdens,threspB=threspB,threstaunorm=threstaunorm ;+ ; Name: FOR_PINPOINT ; ; Purpose: determine weighted center of mass along LOS using WL polarization ratio ; ; INPUTS ; ; POL == polarization ration (pB/B) ; dimension POS ; PB = polarization brightness ; dimension POS ; R3D -- radial distance along LOS ; dimension LOS X POS ; TAU3D -- integration angle from TS along LOS ; dimension LOS X POS ; DENSITY -- density ; dimension LOS X POS ; RPOS -- plane of sky radius (relates to elongation) ; dimension POS ; DISTOBS -- distance to observer (defaults to 215 Rsun) ; ; THRESDENS, THRESPB, THRESTAUNORM being tested-- at the moment are only internal to this program, ; to be set by uncommenting below ; ; OUTPUTS ; XX_PLUS -- X position of solution with mass in front of Thomson sphere ; XX_MINUS -- X position of solution with mass behind Thomson sphere ; XX_GROUND -- X position of ground truth solution - position center of mass along LOS ; TAU_PLUS -- Tau position along LOS of solution with mass in front of Thomson sphere ; TAU_MINUS -- Tau position along LOS of solution with mass behind Thomson sphere ; TAU_GROUND -- Tau position along LOS of ground truth solution - position center of mass along LOS ; ; Called by: FOR_INTEGRATE ; ; Calls: FOR_BFUN, FOR_CFUN ; ; ; HISTORY: ; Written by Sarah Gibson 2020 ; ; December 2021-- made DISTOBS a required input ; January 2022 -- added TPOLF, TPOLB, TPOLG angular distance from TS functionality ; also changed the way center of mass was calculated, making it ; relative to observer not TS ; July 2024 passed through ulimb ; Jan 2025 -- fixed bug where large epsilon (close to pi/2) ; had bad fit -- replaced with point source ; also passed through los and losuse to be able to make the integral ; work for both TLOS and TAU ; April 2025 -- updated integral of TPOLG to use int_tabulated ;- default,dimrad,201 ;default,threstaunorm,.01 ;default,thresdens,5d7 ;default,threspB,1d-6 ; distance to observer is usually 1 AU ; set in distobs ; elongation angle epslon epslon = atan(rpos,distobs) ; points of intersection with Thomson sphere ; this is d in paper notation rp2=rpos*cos(epslon) ; Br/Bt = PR PR=(1-pol)/(1+pol) ; from Billings/Guhathakurta thesis ; Br/Bt = 1. - sin^2(chi)*scatfunc(r) ; chi=asin(sqrt(scatfunc*(1.-PR))) ; ; since we don't know r, we can assume the sun as a point source, ; scatfunc --> 1, this reduces to ; chi=asin(sqrt(1-PR)) ; chi=acos(sqrt(PR)) ; ; this is what is in DeForest et al eq 3 and ok for > 2-3 Rsun ; BUT introduces error below this ; For lower heights ; we will introduce another dimension and solve ; for best fit to radius - knowing that it is at minimum rpos ; and we will assume it lies within 2.*rpos of the plane of sky ; *this intrinsically assumes that we can treat the scattering ; *as a single super-particle at a particular position ; *but then so does this entire analysis ; get=size(rpos) dimrpos=get[1] radius=dblarr(dimrpos,dimrad) PR_use=radius*0. epslon_use=radius*0. for i = 0,(dimrad-1) do begin argument=(rpos[*]^2+((i-(dimrad-1)/2.)*2.*rpos[*]/((dimrad-1)/2.))^2) radius[*,i]=sqrt(argument) PR_use[*,i]=PR epslon_use[*,i]=epslon endfor scatfunc=for_bfun(radius,ulimb=ulimb)/for_cfun(radius,ulimb=ulimb) ; also consider both foreground and background solutions ; be careful of how scattering angle chi is defined -- ; in calculations using chi below I'm assuming ; chi_minus is greater than 90, thus 180-chi_plus ; If one considers the two signs of the sqrt, one gets the two angles defined ; as +/- each other because sin(-x)=-sin(x) ; which is fine, but would be treated differently below chi_plus=asin(sqrt(scatfunc*(1.-PR_use))) chi_minus=!dpi - chi_plus ; ; note -- scatfunc*(1.-PR) may be negative -- this just means ; the superparticle definitely is not located in such points ; note- - there may be cases where 2*rpos is not enough ; and as a result no solution will be found along the entire LOS. ; These will be commented on at end of program (bad sim points) ; ***CHANGED THIS - this tends to happen for large epsilon ; so no need for this fit there (point source is ok) ; ; point source chi_plus_ps=acos(sqrt(PR)) chi_minus_ps= !dpi - chi_plus_ps ; ; from poster ; feature in the foreground ; ; pi/2-xi_plus (which is theta) + (pi-chi_plus) (which is chi_plus') + epslon=pi xi_plus = epslon_use - chi_plus + !dpi/2.d0 xi_plus_ps = epslon - chi_plus_ps + !dpi/2.d0 ; ; feature in background ; (pi/2-xi_minus)(where xi_minus is going to be smaller than xi_plus, or even negative if behind the POS) ; + (pi-chi_minus) (which is chi_minus') + epslon = pi ;xi_minus = -!dpi/2.d0 + chi_minus - epslon_use ;xi_minus_ps = -!dpi/2.d0 + chi_minus_ps - epslon ; note I had a confusing but at least consistently wrong sign convention that ; I am replacing with a more correct but same result convention, here and below xi_minus = epslon_use - chi_minus + !dpi/2.d0 xi_minus_ps = epslon - chi_minus_ps + !dpi/2.d0 ; ; distance from Earth to localized density ; dist_plus=distobs*(sin(!dpi/2.d0 - xi_plus)/sin(!dpi-chi_plus)) dist_plus_ps=distobs*(sin(!dpi/2.d0 - xi_plus_ps)/sin(!dpi-chi_plus_ps)) ;dist_minus=distobs*(sin(!dpi/2.d0 + xi_minus)/sin(!dpi-chi_minus)) ;dist_minus_ps=distobs*(sin(!dpi/2.d0 + xi_minus_ps)/sin(!dpi-chi_minus_ps)) ; sign convention adjusted as above dist_minus=distobs*(sin(!dpi/2.d0 - xi_minus)/sin(!dpi-chi_minus)) dist_minus_ps=distobs*(sin(!dpi/2.d0 - xi_minus_ps)/sin(!dpi-chi_minus_ps)) ; ; these next will be the same ; they are radial distance to scattering points ; rad_plus=distobs*(sin(epslon_use)/sin(!dpi-chi_plus)) rad_plus_ps=distobs*(sin(epslon)/sin(!dpi-chi_plus_ps)) rad_minus=distobs*(sin(epslon_use)/sin(!dpi-chi_minus)) rad_minus_ps=distobs*(sin(epslon)/sin(!dpi-chi_minus_ps)) ; ; now collapse on to the correct radii xi_plus_use=rpos*0. xi_minus_use=rpos*0. rad_plus_use=rpos*0. rad_minus_use=rpos*0. dist_plus_use=rpos*0. dist_minus_use=rpos*0. xx_plus=rpos*0. xx_minus=rpos*0. tau_plus=rpos*0. tau_minus=rpos*0. xx_plus_nocor=rpos*0. xx_minus_nocor=rpos*0. tau_plus_nocor=rpos*0. tau_minus_nocor=rpos*0. xx_ground=rpos*0. tau_ground=rpos*0. dist_ground=rpos*0. dist_fromTS_ground=rpos*0. get=size(r3D) dimlos=get[1] xx3D=dblarr(dimlos) dfromTS=dblarr(dimlos) dist3D=dblarr(dimlos) xi3D=dblarr(dimlos) for i = 0,dimrpos-1 do begin xi3D[*] = epslon[i] + tau3D[*,i] xx3D[*]=r3D[*,i]*cos(!dpi/2.d0 -xi3D) ; ; note, dist3D is now distance from Earth and always positive ; script-l in paper notation ; dfromTS is distance from TS ; s in paper notation ; both are along LOS ; and XX3D is distance from plane of sky along X axis ; dfromTS[*]=r3D[*,i]*sin(tau3D[*,i]) dearthTS=sqrt(distobs^2-rp2[i]^2) dist3D[*]=dearthTS - dfromTS[*] ; for low elongations, parallel LOS, dfromTS and XX are the same ; **NEW** ; only bother to do this if epslon gt 0.1 (point source) ; if epslon[i] lt 0.1 then begin if pol[i]*0. eq 0. and pol[i] ne -8888. then begin posrad=where(abs(rad_plus[i,*]-radius[i,*]) eq min(abs(rad_plus[i,*]-radius[i,*]))) ; ; there should be two solutions -- symmetric about the POS *do I mean TS?* ; it should not matter which is chosen ; for comparison, will also calculate distance if we assumed point source ; xi_plus_use[i]=xi_plus[i,posrad[0]] xi_minus_use[i]=xi_minus[i,posrad[0]] rad_plus_use[i]=rad_plus[i,posrad[0]] rad_minus_use[i]=rad_minus[i,posrad[0]] dist_plus_use[i]=dist_plus[i,posrad[0]] dist_minus_use[i]=dist_minus[i,posrad[0]] endif else begin xx_plus[i]=-8888. xx_minus[i]=-8888. xx_ground[i]=-8888. tau_plus[i]=-8888. tau_minus[i]=-8888. tau_ground[i]=-8888. dist_ground[i]=-8888. dist_fromTS_ground[i]=-8888. endelse endif else begin xi_plus_use[i]=xi_plus_ps[i] xi_minus_use[i]=xi_minus_ps[i] rad_plus_use[i]=rad_plus_ps[i] rad_minus_use[i]=rad_minus_ps[i] dist_plus_use[i]=dist_plus_ps[i] dist_minus_use[i]=dist_minus_ps[i] endelse xx_plus[i]=rad_plus_use[i]*cos(!dpi/2.d0 -xi_plus_use[i]) ; these next will be PS all the way xx_plus_nocor[i]=rad_plus_ps[i]*cos(!dpi/2.d0 -xi_plus_ps[i]) xx_minus[i]=rad_minus_use[i]*cos(!dpi/2.d0 -xi_minus_use[i]) xx_minus_nocor[i]=rad_minus_ps[i]*cos(!dpi/2.d0 -xi_minus_ps[i]) ; xx_minus[i]=-rad_minus_use[i]*cos(!dpi/2.d0 -xi_minus_use[i]) ; xx_minus_nocor[i]=-rad_minus_ps[i]*cos(!dpi/2.d0 -xi_minus_ps[i]) ; adjusted sign convention as above tau_plus[i]= xi_plus_use[i] - epslon[i] tau_minus[i]= xi_minus_use[i] - epslon[i] tau_plus_nocor[i]= xi_plus_ps[i] - epslon[i] tau_minus_nocor[i]= xi_minus_ps[i] - epslon[i] ; ; for ground truth check, figure out where along the line of sight ; the center of mass is for distributed source ; relative to intersection with TS ; ; allow exclusion of density below a given threshold ; denslos=Density[*,i] taulos=tau3d[*,i] if keyword_set(thresdens) then begin thresground=where(Density[*,i] lt thresdens) if min(thresground) ne -1 then denslos[thresground]=0. endif ; I tried to do COM integrals over tau but ran into problems ; - safer to do over distance but be careful about dlos ; which will not be constant for TAU ; if strupcase(losuse) eq 'TAU' then begin ; dist_ground[i]=total(denslos*dist3D*rp2[i]/cos(taulos)^2)/total(denslos*rp2[i]/cos(taulos)^2) ; below is same just rewritten differently ; dist_fromTS_ground[i]=rp2[i]*total(denslos*sin(taulos)/cos(taulos)^3)/total(denslos/cos(taulos)^2) integrand1=denslos*dfromts*rp2[i]/cos(taulos)/cos(taulos) ; integrand1=denslos*rp2[i]*rp2[i]*tan(taulos)/cos(taulos)/cos(taulos) integrand2=denslos*rp2[i]/cos(taulos)/cos(taulos) losp=taulos endif else begin ; dist3D is script-l ; dfromTS is s (also same as dearthTS - dist3D) ; rp2 is d in paper ; TLOS - constant dlos so it cancels outside the integrand ; dist_ground[i]=total(denslos*dist3D)/total(denslos) ; this next is Scom ; dist_fromTS_ground[i]=dearthTS-dist_ground[i] ; integrand1=denslos*dfromts integrand2=denslos losp=dfromTS endelse ix=uniq(reform(losp)) if (n_elements(ix) gt 1) then begin dist_fromTS_ground[i]=int_tabulated(reform(losp[ix]),reform(integrand1[ix])) dist_fromTS_ground[i]=dist_fromTS_ground[i]/int_tabulated(reform(losp[ix]),reform(integrand2[ix])) endif else print,'badlos for center of mass' tau_ground[i]=atan(dist_fromTS_ground[i],rp2[i]) r_ground=sqrt(dist_fromTS_ground[i]^2+rp2[i]*rp2[i]) xx_ground[i]=r_ground*sin(epslon[i]+tau_ground[i]) ; ; if abs(xx_ground[i]-total(denslos*xx3D)/total(denslos)) gt 1d-8 then stop ; - so this could have been shortened to ; xx_ground[i]=total(denslos*xx3D)/total(denslos) ; for background, non parallel lines of sight ; means that far field will be intersecting huge heights ; along entire LOS and very little (if any) points behind POS ; and low densities if so, so XPOLG will be positive (and big) ; which is why TPOLG is probably more useful ; ; print,'ground truth center of mass distance from plane of sky=',xx_ground[i] ; print,'test ground truth=',total(denslos*xx3D)/total(denslos) ; print,'min/max distance from plane of sky - ground truth',minmax(xx3D) ; print,'polarization position (foreground)=',xx_plus[i] ; print,'polarization position (foreground) - point source assumption =',xx_plus_nocor[i] ; print,'polarization position (background)=',xx_minus[i] ; print,'polarization position (background) - point source assumption =',xx_minus_nocor[i] endfor ; ; try removing background by putting a threshold on pB ; if keyword_set(threspB) then begin threspol=where(pB lt threspB) if min(threspol) ne -1 then begin xx_plus[threspol]=-8888. xx_minus[threspol]=-8888. tau_plus[threspol]=-8888. tau_minus[threspol]=-8888. endif endif if keyword_set(threstaunorm) then begin threstau_plus=where(tau_plus lt threstaunorm*max(abs(tau_plus))) if min(threstau_plus) ne -1 then begin xx_plus[threstau_plus]=-8888. tau_plus[threstau_plus]=-8888. endif threstau_minus=where(tau_minus lt threstaunorm*max(abs(tau_minus))) if min(threstau_minus) ne -1 then begin xx_minus[threstau_minus]=-8888. tau_minus[threstau_minus]=-8888. endif endif ; ; there can be places where we get NaNs because ; we haven't extended radius far enough. For now just ; report on these and make them bad points ; **should have removed these by using PS for higher elongation test=where(xx_plus*0. ne 0. or xx_minus*0. ne 0.) if min(test) ne -1 then begin print,n_elements(test),' bad sim points' ; print,xx_plus[test] ; print,xx_minus[test] xx_plus[test]=-8888. xx_minus[test]=-8888. tau_plus[test]=-8888. tau_minus[test]=-8888. endif test=where(xx_ground*0. ne 0.) if min(test) ne -1 then begin print,n_elements(test),' bad ground points' ; print,xx_ground[test] xx_ground[test]=-8888. tau_ground[test]=-8888. endif ; ; uncomment ; if want to check point source vs numerically solved ; ;tau_plus=tau_plus_nocor ;tau_minus=tau_minus_nocor ;xx_plus=xx_plus_nocor ;xx_minus=xx_minus_nocor end