wcsxy2sph.pro

basic median filter simulation

    temp = where(larger ne 0.d0, Ntemp)
    if Ntemp GT 0 then ratio[temp] = sin(pi/1.2d1*smaller[temp]/larger[temp])/$
                      (cos(pi/1.2d1*smaller[temp]/larger[temp]) - sqrt(0.5d0))

    temp = where(larger eq 0.d0, Ntemp)
    if Ntemp GT 0 then ratio[temp] = 1.d0
    rho = 1.d0 - (larger)^2*(1.d0 - 1.d0/sqrt(2.d0 + ratio^2))

    temp = where((abs(xx) gt abs(yy)) and (ratio ne 0.d0), Ntemp)
    if Ntemp GT 0 then ratio[temp] = 1.d0/ratio[temp]

    temp = where((abs(xx) gt abs(yy)) and (ratio eq 0.d0), Ntemp)
; use a kludge to produce the correct value for 1/0 without generating an error 
    if Ntemp GT 0 then ratio[temp] = tan(pi2)

    for i = 0l, n_x-1 do begin
      case face[i] of
        0:begin
          if (xx[i] ne 0.d0) then phi[i] = atan(-ratio[i]) else $
          if (yy[i] le 0.d0) then phi[i] = 0.d0 else $
          if (yy[i] gt 0.d0) then phi[i] = pi

          if (yy[i] ne 0.d0) then theta[i] = asin(rho[i]) else $
          if (xx[i] le 0.d0) then theta[i] = -pi2 else $
          if (xx[i] gt 0.d0) then theta[i] = pi2

          if (yy[i] gt 0.d0) then begin
            if (xx[i] lt 0.d0) then phi[i] = phi[i] - pi $
            else if (xx[i] gt 0.d0) then phi[i] = phi[i] + pi
          endif
        end
        1:begin
          if (xx[i] ne 0.d0) then begin
            if (yy[i] ne 0.d0) then $
             phi[i] = xx[i]/abs(xx[i])*acos(sqrt(rho[i]^2*(1.d0 + ratio[i]^2)/$
                             (ratio[i]^2 + rho[i]^2))) $
            else phi[i] = xx[i]/abs(xx[i])*acos(rho[i]) 
          endif else phi[i] = 0.d0
          if (yy[i] ne 0.d0) then theta[i] = yy[i]/abs(yy[i])*acos(rho[i]/$
                                            cos(phi[i])) else theta[i] = 0.d0
        end
        2:begin
          if (yy[i] ne 0.d0) then begin
            if (xx[i] gt 0.d0) then $
               phi[i] = pi - asin(sqrt(rho[i]^2*(1.d0 + ratio[i]^2)/$
                             (ratio[i]^2 + rho[i]^2))) $
            else if (xx[i] lt 0.d0) then $
               phi[i] = asin(sqrt(rho[i]^2*(1.d0 + ratio[i]^2)/$
                             (ratio[i]^2 + rho[i]^2))) $
            else phi[i] = pi2
            theta[i] = yy[i]/abs(yy[i])*acos(rho[i]/abs(sin(phi[i])))
          endif else begin 
            theta[i] = 0.d0
            if (xx[i] gt 0.d0) then phi[i] = pi - asin(rho[i]) $
            else if (xx[i] lt 0.d0) then phi[i] = asin(rho[i]) $
            else phi[i] = pi2
          endelse
        end
        3:begin
          if (yy[i] ne 0.d0) then begin
            if (xx[i] gt 0.d0) then $
              phi[i] = acos(sqrt(rho[i]^2*(1.d0 + ratio[i]^2)/$
                            (ratio[i]^2 + rho[i]^2))) - pi $
            else if (xx[i] lt 0.d0) then $
              phi[i] = -acos(sqrt(rho[i]^2*(1.d0 + ratio[i]^2)/$
                       (ratio[i]^2 + rho[i]^2))) + pi $
            else phi[i] = pi
            theta[i] = yy[i]/abs(yy[i])*acos(-rho[i]/cos(phi[i])) 
          endif else begin
            theta[i] = 0.d0
            if (xx[i] gt 0.d0) then phi[i] = acos(rho[i]) - pi $
            else if (xx[i] lt 0.d0) then phi[i] = -acos(rho[i]) + pi $
            else phi[i] = pi
          endelse
        end
        4:begin
          if (yy[i] ne 0.d0) then begin
            if (xx[i] gt 0.d0) then $
               phi[i] = -asin(sqrt(rho[i]^2*(1.d0 + ratio[i]^2)/$
                             (ratio[i]^2 + rho[i]^2))) $
            else if (xx[i] lt 0.d0) then $
               phi[i] = asin(sqrt(rho[i]^2*(1.d0 + ratio[i]^2)/$
                             (ratio[i]^2 + rho[i]^2))) - pi $
            else phi[i] = -pi2
            theta[i] = yy[i]/abs(yy[i])*acos(-rho[i]/sin(phi[i]))
          endif else begin
            theta[i] = 0.d0
            if (xx[i] gt 0.d0) then phi[i] = -asin(rho[i]) $
            else if (xx[i] lt 0.d0) then phi[i] = asin(rho[i]) - pi $
            else phi[i] = -pi2
          endelse
        end
        5:begin
          if (xx[i] ne 0.d0) then phi[i] = atan(ratio[i]) $
          else if (yy[i] le 0.d0) then phi[i] = pi $
          else if (yy[i] gt 0.d0) then phi[i] = 0.d0

          if (yy[i] ne 0.d0) then theta[i] = asin(-rho[i]) $
          else if (xx[i] le 0.d0) then theta[i] = -pi2 $
          else if (xx[i] gt 0.d0) then theta[i] = pi2

          if (yy[i] lt 0.d0) then begin
            if (xx[i] lt 0.d0) then phi[i] = phi[i] - pi $
            else if (xx[i] gt 0.d0) then phi[i] = phi[i] + pi
          endif
        end
      endcase
    endfor
  end
  
  'TSC':begin

    xx=xx/45.0d0
    yy=yy/45.0d0
;
;   If the faces are not defined, assume that the faces need to be defined
;   and the whole sky is displayed as a "sideways T".
;
        if noface eq 1 then begin

                face=intarr(n_elements(xx))

                face1 = where((xx le 1.0) and (yy le 1.0) and (yy ge -1.0),nf1)
                if nf1 gt 0 then begin
                        face[face1]=1
                endif

                face4 = where((xx gt 5.0),nf4)
                if nf4 gt 0 then begin
                        face[face4]=4
                        xx[face4]=xx[face4]-6.0d0
                endif

                face3 = where((xx le 5.0) and (xx gt 3.0),nf3)
                if nf3 gt 0 then begin
                        face[face3]=3
                        xx[face3]=xx[face3]-4.0d0
                endif

                face2 = where((xx le 3.0) and (xx gt 1.0),nf2)
                if nf2 gt 0 then begin
                        face[face2]=2
                        xx[face2]=xx[face2]-2.0d0
                endif

                face0 = where((xx le 1.0) and (yy gt 1.0),nf0)
                if nf0 gt 0 then begin
                        face[face0]=0
                        yy[face0]=yy[face0] - 2.0
                endif

                face5 = where((xx le 1.0) and (yy lt -1.0),nf5)
                if nf5 gt 0 then begin
                        face[face5]=5
                        yy[face5]=yy[face5] + 2.0
                endif

        endif

    rho = sin(atan(1.0d0/sqrt(xx^2 + yy^2)))
    phi = double(x - x)
    theta = double(x - x)
    for i = 0l, n_x - 1 do begin
      case face[i] of
        0:begin
          phi[i] = atan(xx[i],-yy[i])
          theta[i] = asin(rho[i])
        end
        1:begin
          if (xx[i] ne 0.d0) then begin
            if (xx[i] ge 0.d0) then $
             phi[i] = atan(sqrt((1.d0/rho[i]^2- 1.d0)/(1 + (yy[i]/xx[i])^2))) $
            else phi[i] =atan(-sqrt((1.d0/rho[i]^2 - 1.d0)/$
                               (1 + (yy[i]/xx[i])^2)))
            theta[i] = atan(yy[i]/xx[i]*sin(phi[i]))
          endif else begin
            phi[i] = 0.d0
            if (yy[i] ge 0.d0) then theta[i] = acos(rho[i]) $
            else theta[i] = -acos(rho[i])
          endelse
        end
        2:begin
; The point theta = 0, phi = Pi/2 lies in this region, allowing 
; rho = Cos[theta]*Sin[phi] to be 1, causing an infinite quantity in the
; equation for phi
          if (rho[i] eq 1.d0) then begin
            phi[i] = pi2
            theta[i] = 0.d0
          endif else if (xx[i] gt 1.d-14) then begin
           phi[i] = atan(-sqrt((1.d0 + (yy[i]/xx[i])^2)/$
                                (1.d0/rho[i]^2 - 1.d0)))+pi
            theta[i] = atan(-yy[i]/xx[i]*cos(phi[i]))
          endif else if (xx[i] lt -1.d-14) then begin
            phi[i]=atan(sqrt((1.d0+(yy[i]/xx[i])^2)/(1.d0/rho[i]^2 - 1.d0)))
            theta[i] = atan(-yy[i]/xx[i]*cos(phi[i]))
          endif else begin
             phi[i] = pi2
            if (yy[i] ge 0) then theta[i] = acos(rho[i]/sin(phi[i])) $
            else theta[i] = -acos(rho[i]/sin(phi[i]))
          endelse
        end
        3:begin
          if (abs(xx[i]) ge 1.d-5) then begin
            if (xx[i] gt 0.d0) then $
           phi[i] = atan(sqrt((1.d0/rho[i]^2 - 1.d0)/$
                          (1 + (yy[i]/xx[i])^2)))-pi $
        else phi[i] = atan(-sqrt((1.d0/rho[i]^2 - 1.d0)/$
                            (1 + (yy[i]/xx[i])^2)))+pi
            theta[i] = atan(-yy[i]/xx[i]*sin(phi[i]))
          endif else begin
            if (xx[i] ge 0.d0) then phi[i] = -pi $
            else phi[i] = pi
            if (yy[i] ge 0) then theta[i] = acos(rho[i]) $
            else theta[i] = -acos(rho[i])
          endelse
        end
        4:begin
          if (rho[i] eq 1.d0) then begin
            phi[i] = -pi2
            theta[i] = atan(yy[i]/xx[i])
          endif else if (xx[i] gt 1.d-14) then begin
           phi[i]=atan(-sqrt((1.d0 + (yy[i]/xx[i])^2)/(1.d0/rho[i]^2 - 1.d0)))
            theta[i] = atan(yy[i]/xx[i]*cos(phi[i]))
          endif else if (xx[i] lt -1.d-14) then begin
            phi[i]=atan(sqrt((1.d0+(yy[i]/xx[i])^2)/(1.d0/rho[i]^2 - 1.d0)))-pi
            theta[i] = atan(yy[i]/xx[i]*cos(phi[i]))
          endif else begin
             phi[i] = 1.5d0*!pi
            if (yy[i] ge 0) then theta[i] = acos(rho[i]) $
            else theta[i] = -acos(rho[i])
          endelse
        end
        5:begin
          phi[i] = atan(xx[i],yy[i])
          theta[i] = asin(-rho[i])
        end

      endcase
    endfor
  end
  
  'HPX':begin
    hpx_k = 3.D ; The main HEALPIX parameters (see Calabretta 2007, MNRAS)
    hpx_h = 4.D ;
    ylim = 90D *(hpx_k-1)/hpx_h
    phi=xx*1.D
    theta=yy*1.D
    
    eqfaces = where( abs(yy) le ylim, complement=polfaces)

    ; equatorial region
    if eqfaces[0] ne -1 then begin
        phi[eqfaces]=xx[eqfaces]/radeg
        theta[eqfaces]=asin(yy[eqfaces]*hpx_h/90.D/hpx_k)
	hpx_bad = where(xx[eqfaces] lt -180 or xx[eqfaces] gt 180)
	if hpx_bad[0] ne -1 then begin
		phi[eqfaces[hpx_bad]]=!VALUES.F_NAN
		theta[eqfaces[hpx_bad]]=!VALUES.F_NAN
	end

    endif

    ;polar regions
    if polfaces[0] ne -1 then begin
        hpx_sig = (hpx_k+1)/2.D -abs(yy[polfaces]*hpx_h)/180.D
        hpx_omega = ((hpx_k mod 2 eq 1) or yy[polfaces] gt 0)*1.D
        hpx_xc=-180+(2*floor((xx[polfaces]+180)*hpx_h/360+(1-hpx_omega)/2)+hpx_omega)*180.D/hpx_h
        phi[polfaces] = (hpx_xc + (xx[polfaces]-hpx_xc)/hpx_sig)/radeg
        theta[polfaces] = ((yy[polfaces] gt 0)*2-1)*asin(1-hpx_sig^2/hpx_k)

	; points outside reasonable area
	hpx_bad = where( ((abs(((xx[polfaces]-hpx_xc)*hpx_h) mod 180)+(abs(yy[polfaces])-ylim)*hpx_h) gt 45*(3+hpx_h-hpx_k)) or (xx[polfaces] lt -180) or (xx[polfaces] gt 180))
	if hpx_bad[0] ne -1 then begin
		phi[polfaces[hpx_bad]]=!VALUES.F_NAN
		theta[polfaces[hpx_bad]]=!VALUES.F_NAN
	end
    endif
  end
  else:message,strupcase(projection_type) + $
               ' is not a valid projection type.  Reset CTYPE1 and CTYPE2'

endcase

; Convert from "native" coordinate system to "standard" coordinate system
; if the CRVAL keyword is set.  Otherwise, assume the map projection is 
; complete 

 phi = phi*radeg
 theta = theta*radeg
 
 if ( N_elements(crval) GE 2 ) then begin

 

  if (n_elements(longpole) eq 0) then longpole = 1.8d2

  if N_elements(map_type) EQ 0 then $
           map_type = where(projection_type EQ map_types)
   map_type = map_type[0]
   conic = (map_type GE 13) and (map_type LE 16) 
   zenithal =  ((map_type GE 1) and (map_type LE 8)) or (map_type EQ 26) 
   if conic then theta0 = pv2_1 else if zenithal then theta0 = 90 $
            else theta0 = 0 
   wcs_rotate, longitude, latitude, phi, theta, crval, longpole=longpole, $
           theta0 = theta0, latpole = latpole, /REVERSE
 endif else begin    ;no rotation from standard to native coordinates

  latitude = theta
  longitude = phi

endelse

; CONVERT LONGITUDE FROM -180 TO 180 TO 0 TO 360

 temp = where(longitude lt 0.d0, Nneg)
 if (Nneg GT 0) then longitude[temp] = longitude[temp] + 3.6d2
 temp = where(longitude ge 3.6d2-1.d-2, Nneg)
 if (Nneg GT 0) then longitude[temp] = longitude[temp] - 3.6d2


 return
 end