wcsxy2sph.pro

basic median filter simulation

    if (pv2_1 eq 0) then message,$
      'pv2_1 = 0 for BON map projection is better done with SFL map projection'
    theta_1 = pv2_1/radeg
    y_0 = 1.d0/tan(theta_1) + theta_1
    s = theta_1/abs(theta_1)
    theta = y_0 - s*sqrt(xx^2 + (y_0*radeg - yy)^2)/radeg 
    phi = s*(y_0 - theta)*atan(s*xx/(y_0*radeg - theta),$
                               (y_0*radeg - yy)/(y_0*radeg - theta))/cos(theta)
  end
  
  'PCO':begin
; Determine where y = 0 and assign theta to 0 for these points.  The reason
; for doing this separately is that the intial condition for theta in the
; numerical solution is sign(y)*45 which only works for y not = 0.
    bad = where(yy eq 0)
    good = where(yy ne 0)
    theta = double(xx - xx)
    if (bad[0] ne -1) then theta[bad] = 0.d0

; Find theta numerically.
    tolerance = 1.d-11
    tolerance_2 = 1.d-11
    if (good[0] ne -1) then begin
      theta_p = double(xx - xx)
      theta_p[good] = pi2*yy[good]/abs(yy[good])
      theta_n = double(xx - xx)
      f_p = double(xx - xx)
      f_p[good] = xx[good]^2 - 2.d0*radeg*(yy[good] - radeg*theta_p[good])/$
                  tan(theta_p[good]) + (yy[good] - radeg*theta_p[good])^2
      f_n = double(xx - xx) - 999.d0
      lambda = double(xx - xx)
      f = double(xx - xx)
      repeat begin
        case_1 = where((yy ne 0.d0) and (f_n lt (-1.d2)))
        case_2 = where((yy ne 0.d0) and (f_n ge (-1.d2)))
        if (case_1[0] ne -1) then lambda[case_1] = 0.5d0
        if (case_2[0] ne -1) then $
          lambda[case_2] = f_p[case_2]/(f_p[case_2] - f_n[case_2])
        lambda[good] = 1.d-1 > (9.d-1 < lambda[good])
        theta[good] = (1.d0 - lambda[good])*theta_p[good] + $
                      lambda[good]*theta_n[good]
        f[good] = xx[good]^2 - 2.d0*radeg*(yy[good] - radeg*theta[good])/$
                  tan(theta[good]) + (yy[good] - radeg*theta[good])^2
        neg = where((yy ne 0.d0) and (f lt 0.d0))
        pos = where((yy ne 0.d0) and (f gt 0.d0))
        if (neg[0] ne -1) then begin
          f_n[neg] = f[neg]
          theta_n[neg] = theta[neg]
        end
        if (pos[0] ne -1) then begin
          f_p[pos] = f[pos]
          theta_p[pos] = theta[pos]
        end
      endrep until ((max(abs(theta_p - theta_n)) lt tolerance) or $
                    (max(abs(f)) lt tolerance_2))
    endif

; Determine phi differently depending on whether theta = 0 or not.
    bad = where(theta eq 0.d0)
    good = where(theta ne 0.d0)
    phi = double(x - x)
    if (bad[0] ne -1) then phi[bad] = xx[bad]/radeg
    phi[good] = atan(xx[good]/radeg*tan(theta[good]),$
       1.d0 - (yy[good]/radeg - theta[good])*tan(theta[good]))/sin(theta[good])
  end
  
  'SFL':begin
    phi = xx/(radeg*cos(yy/radeg))
    theta = yy/radeg
  end
  
  'PAR':begin
  
    theta = 3.d0*asin(yy/pi/radeg)
    phi = xx/(1.d0 - 4.d0*(yy/pi/radeg)^2)/radeg
  end
  
  'AIT':begin
  z2 = 1.d0 - (xx/(4.d0*radeg))^2 - (yy/(2.d0*radeg))^2
  bad = where(z2 lt 0.5d0,nbad)
  z = sqrt(z2)
  phi = 2.d0*atan(z*xx/(2.d0*radeg),2.d0*z^2 - 1.d0)
  theta = asin(yy*z/radeg)
  if nbad gt 0 then begin
       phi[bad] = !values.d_nan
      theta[bad] = !values.d_nan
   endif

  end
  
  'MOL':begin
    phi = pi*xx/(radeg*2.d0*sqrt(2.d0 - (yy/radeg)^2))
    arg = 2.d0*asin(yy/(sqrt(2.d0)*radeg))/pi + $
                 yy*sqrt(2.d0 - (yy/radeg)^2)/1.8d2
    
   theta = asin(2.d0*asin(yy/(sqrt(2.d0)*radeg))/pi + $
                 yy*sqrt(2.d0 - (yy/radeg)^2)/1.8d2)
    
  end
  
  'CSC':begin
    xx = xx/4.5d1
    yy = yy/4.5d1

;
;   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

; Define array of numerical constants used in determining alpha and beta1.
    p = dblarr(7,7)
    p[0,0] = -0.27292696d0
    p[1,0] = -0.07629969d0
    p[0,1] = -0.02819452d0
    p[2,0] = -0.22797056d0
    p[1,1] = -0.01471565d0
    p[0,2] = 0.27058160d0
    p[3,0] = 0.54852384d0
    p[2,1] = 0.48051509d0
    p[1,2] = -0.56800938d0
    p[0,3] = -0.60441560d0
    p[4,0] = -0.62930065d0
    p[3,1] = -1.74114454d0
    p[2,2] = 0.30803317d0
    p[1,3] = 1.50880086d0
    p[0,4] = 0.93412077d0
    p[5,0] = 0.25795794d0
    p[4,1] = 1.71547508d0
    p[3,2] = 0.98938102d0
    p[2,3] = -0.93678576d0
    p[1,4] = -1.41601920d0
    p[0,5] = -0.63915306d0
    p[6,0] = 0.02584375d0
    p[5,1] = -0.53022337d0
    p[4,2] = -0.83180469d0
    p[3,3] = 0.08693841d0
    p[2,4] = 0.33887446d0
    p[1,5] = 0.52032238d0
    p[0,6] = 0.14381585d0

; Calculate alpha and beta1 using numerical constants
    sum = double(x - x)
    for j = 0,6 do for i = 0,6 - j do sum = sum + p[i,j]*xx^(2*i)*yy^(2*j)
    alpha = xx + xx*(1 - xx^2)*sum

    sum = double(x - x)
    for j = 0,6 do for i = 0,6 - j do sum = sum + p[i,j]*yy^(2*i)*xx^(2*j)
    beta1 = yy + yy*(1 - yy^2)*sum

; Calculate theta and phi from alpha and beta1; the method depends on which
; face the point lies on
    phi = double(x - x)
    theta = double(x - x)
    for i = 0l, n_x - 1 do begin
      case face[i] of
        0:begin
          if (beta1[i] eq 0.d0) then begin
            if (alpha[i] eq 0.d0) then begin
              theta[i] = pi2
; uh-oh lost information if this happens
              phi[i] = 0.d0
            endif else begin
              phi[i] = alpha[i]/abs(alpha[i])*pi2
              theta[i] = atan(abs(1.d0/alpha[i]))
            endelse
          endif else begin
            phi[i] = atan(alpha[i],-beta1[i])
            theta[i] = atan(-cos(phi[i])/beta1[i])
          endelse
; ensure that the latitudes are positive
          theta[i] = abs(theta[i])
        end
        1:begin
          phi[i] = atan(alpha[i])
          theta[i] = atan(beta1[i]*cos(phi[i]))
        end
        2:begin
          if (alpha[i] eq 0.d0) then phi[i] = pi2 else $
            phi[i] = atan(-1.d0/alpha[i])
          if (phi[i] lt 0.d0) then phi[i] = phi[i] + pi
          theta[i] = atan(beta1[i]*sin(phi[i]))
        end
        3:begin
          phi[i] = atan(alpha[i])
          if (phi[i] gt 0.d0) then phi[i] = phi[i] - pi else $
          if (phi[i] lt 0.d0) then phi[i] = phi[i] + pi
          theta[i] = atan(-beta1[i]*cos(phi[i]))
        end
        4:begin
          if (alpha[i] eq 0.d0) then phi[i] = -pi2 else $
            phi[i] = atan(-1.d0/alpha[i])
          if (phi[i] gt 0.d0) then phi[i] = phi[i] - pi
          theta[i] = atan(-beta1[i]*sin(phi[i]))
        end
        5:begin
          if (beta1[i] eq 0.d0) then begin
            if (alpha[i] eq 0.d0) then begin
              theta[i] = -pi2
; uh-oh lost information if this happens
              phi[i] = 0.d0
            endif else begin
              phi[i] = -alpha[i]/abs(alpha[i])*pi2 
              theta[i] = -atan(abs(1.d0/alpha[i]))
            endelse
          endif else begin
            phi[i] = atan(alpha[i],beta1[i])
            theta[i] = atan(-cos(phi[i])/beta1[i])
          endelse
; ensure that the latitudes are negative
          theta[i] = -abs(theta[i])
        end

      endcase
    endfor
  end
  
  'QSC':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


; First determine the quadrant in which each points lies.  Calculate the
; ratio (alpha/beta1) for each point depending on the quadrant.  Finally,
; use this information and the face on which the point lies to calculate
; phi and theta.
    theta = double(x - x)
    phi = double(x - x)
    rho = double(x - x)
    ratio = double(x - x)
    larger = double(x - x)
    smaller = double(x - x)

    temp = where(abs(yy) ge abs(xx), Ntemp)
    if Ntemp GT 0 then larger[temp] = yy[temp]
    temp = where(abs(xx) gt abs(yy), Ntemp )
    if Ntemp GT 0 then larger[temp] = xx[temp]
    temp = where(abs(yy) lt abs(xx), Ntemp )
    if Ntemp GT 0 then smaller[temp] = yy[temp]
    temp = where(abs(xx) le abs(yy), Ntemp)
    if Ntemp GT 0 then smaller[temp] = xx[temp]