function radproject,lat,lon,inrad=inrad,maxrad=maxrad,rrng=rrng,rdel=rdel,$ opaq=opaq,arcone=arcone,theta0=theta0,phi0=phi0,$ xout=xout,yout=yout,zout=zout,rproj=rproj,rout=rout,volout=volout,$ verbose=verbose, _extra=e ;+ ;function radproject ; Projects a ray originating from an arbitrary point on an arbitrarily ; oriented sphere onto plane at \infty. Returns a structure containing ; the radial coordinate, and the corresponding x-coords, y-,z-offsets, ; the magnitude of the projection, the corresponding volume of the ; pillbox, as well as all the input values. ; ;syntax ; rstr=radproject(lat,lon,inrad=inrad,maxrad=maxrad,rrng=rrng,rdel=rdel,$ ; opaq=opaq,arcone=arcone,theta0=theta0,phi0=phi0,$ ; xout=xout,yout=yout,zout=zout,rproj=rproj,rout=rout,volout=volout,$ ; verbose=verbose) ; ;parameters ; lat [INPUT; required] latitude on sphere at which ray is anchored ; lon [INPUT; required] longitude on sphere at which ray is anchored ; * LAT,LON must be scalar ; * (LAT,LON)=(0,0) projects to (Y,Z)=(0,0) ; ;keywords ; inrad [INPUT; default=1] minimum r value of ray to consider ; maxrad [INPUT; default=10] maximum r value of ray to consider ; rrng [INPUT; default=[INRAD,MAXRAD]] filter the output by limiting ; the radii to consider ; * minmax(RRNG) is forced to not spill beyond [INRAD,MAXRAD] ; * use this only if you want everything outside this range to be ; explicitly thrown away. discarded. not zeroed, but removed ; from the output. ; rdel [INPUT; default=0.1] steps in r ; opaq [INPUT; default=0] decides how to treat those portions of rays ; that are occluded by the sphere defined by radius INRAD ; - by default, this is assumed to be completely transparent ; - set OPAQ to negative number to set the opacity as exp(-OPAQ) ; - any value OPAQ>0 results in a completely opaque inner sphere ; (so set /OPAQ to get an opaque inner sphere) ; arcone [INPUT; default=1] area of the opening cone at radius=1 [deg^2] ; theta0 [INPUT; default=0] tilt of sphere [deg] ; phi0 [INPUT; default=0] twist of sphere [deg] ; * theta0 and phi0 are defined in a coordinate frame where the ; projection plane is at x=\infty. ; * it is assumed that to get to the inclined sphere from the ; observer's frame, you first apply THETA0 and then PHI0. ; rout [OUTPUT] radii at which projections are computed ; xout [OUTPUT] computed X ; yout [OUTPUT] projected Y ; zout [OUTPUT] projected Z ; rproj [OUTPUT] projection of the radial unit vector towards the Y-Z ; plane, weighted by the opacity, at each ROUT ; volout [OUTPUT] volume of the pillbox at each ROUT ; verbose [INPUT; default=0] controls chatter ; _extra [JUNK] here only to prevent crashing the program ; ;example ; .run radproject ; ;history ; vinay kashyap (2012dec; supersedes pathrofile) ; added call to SPHTRIG_TRANSLB (VK; 2013feb) ; bug fixes (theta0,phi0=0) (VK; 2013mar) ;- ; usage ok='ok' & np=n_params() & nl=n_elements(lat) & nb=n_elements(lon) if np lt 2 then ok='Insufficient parameters' else $ if nl eq 0 then ok='LAT is not defined' else $ if nb eq 0 then ok='LON is not defined' else $ if nl gt 1 then ok='LAT must be a scalar' else $ if nb gt 1 then ok='LON must be a scalar' if ok ne 'ok' then begin print,'Usage: rstr=radproject(lat,lon,inrad=inrad,maxrad=maxrad,rrng=rrng,rdel=rdel,$' print,' opaq=opaq,arcone=arcone,theta0=theta0,phi0=phi0,$' print,' xout=xout,yout=yout,zout=zout,rproj=rproj,rout=rout,$ print,' verbose=verbose)' print,' projects ray from sphere onto Y-Z plane at \infty' if np ne 0 then message,ok,/informational return,-1L endif ; figure out input keywords vv=0L & if keyword_set(verbose) then vv=long(verbose[0])>1L r0=1.D & if keyword_set(inrad) then r0=abs(double(inrad[0])) r1=10.D & if keyword_set(maxrad) then r1=abs(double(maxrad[0])) radrng=[r0,r1] & nrrng=n_elements(rrng) if nrrng gt 0 then begin if nrrng eq 1 then begin if abs(rrng[0]-r0) lt abs(r1-rrng[0]) then $ radrng[0]=rrng[0]>r0 else radrng[1]=rrng[0]r0 radrng[1]=rr1