wcssph2xy.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
    s = theta_1/abs(theta_1)
    y_0 = 1.d0/tan(theta_1) + theta_1
    a = phi*cos(theta)/(y_0 - theta)
    x = radeg*(y_0 - theta)*sin(a)
    y = radeg*(y_0 - (y_0 - theta)*cos(a))
  end
  
  'PCO':begin
; The equations for x and y are poorly behaved for theta = 0.  Avoid this by
; explicitly assigning values for x and y when theta = 0.
    zero_ind = where(theta eq 0, Nzero)

; create x and y with same structure as longitude
    x = lng*0  & y = x
    if (Nzero GT 0) then begin
      x[zero_ind] = radeg*phi[zero_ind]
      y[zero_ind] = 0.d0
    endif
    good_ind = where(theta ne 0, Ngood)
    if Ngood GT 0 then begin
    x[good_ind] = radeg*sin(phi[good_ind]*sin(theta[good_ind]))/$
                  tan(theta[good_ind])
    y[good_ind] = radeg*(theta[good_ind]+$
        (1.d0 - cos(phi[good_ind]*sin(theta[good_ind])))/tan(theta[good_ind]))
    endif
  end
  
  'SFL':begin
    x = radeg*phi*cos(theta)
    y = radeg*theta
  end
  
  'PAR':begin
    x = radeg*phi*(2.d0*cos(2.d0*theta/3.d0) - 1.d0)
    y = 180.0*sin(theta/3.d0)
  end
  
  'AIT':begin
    alpha = radeg*sqrt(2.d0/(1.d0 + cos(theta)*cos(0.5d0*phi)))
    x = 2.d0*alpha*cos(theta)*sin(0.5d0*phi)
    y = alpha*sin(theta)
  end
  
  'MOL':begin
; Use Newton's method to find a numerical solution to the equation:
;  alpha + 1/2*sin(2*alpha) - 1/2*pi*sin(theta) = 0
    tolerance = 1.0d-14
    alpha = lng*0
    repeat begin
    alpha_old = alpha
    alpha = alpha_old - (alpha_old + 0.5*sin(2.d0*alpha_old) - $
            0.5*pi*sin(theta))/(1.d0 + cos(2.d0*alpha_old))
    endrep until (max(abs(alpha - alpha_old)) lt tolerance)

    x = 2.d0^1.5*phi*radeg*cos(alpha)/pi
    y = sqrt(2.d0)*radeg*sin(alpha)
  end
  
  'CSC':begin
; calculate direction cosines
    l = cos(theta)*sin(phi)
    m = cos(theta)*cos(phi)
    n = sin(theta)

; determine the face on which the x and y coordinates will reside by setting
; rho equal to the maximum of n,m,l,-m,-l,-n which corresponds to faces 0
; through 5 respectively
    rho =  lng*0
    if size(lng,/N_dimen) EQ 0 then  face = 0 else face = lonarr(n_long)

; use an array to store a remapping of the direction cosines.  This way, faces
; 0 and 5 take points on their borders with faces 1-4.  The reason for this is
; that if the max function sees identical values in an array, it takes the 
; index of the first occurrence of that value.
    remap = [0,5,2,1,4,3]

    for i = 0l, n_long-1 do begin
      dir_cos = float([n[i],-n[i],l[i],m[i],-l[i],-m[i]])
      rho[i] = max(dir_cos,temp)
      face[i] = remap[temp]
    endfor

; based on the face determined for each point, find the parameters alpha and
; beta1
    alpha = lng*0
    beta1 = alpha
    for i = 0l, n_long-1 do begin
      case face[i] of
        0:begin
          alpha[i] = l[i]/n[i]
          beta1[i] = -m[i]/n[i]
        end
        1:begin
          alpha[i] = l[i]/m[i]
          beta1[i] = n[i]/m[i]
        end
        2:begin
          alpha[i] = -m[i]/l[i]
          beta1[i] = n[i]/l[i]
        end
        3:begin
          alpha[i] = l[i]/m[i]
          beta1[i] = -n[i]/m[i]
        end
        4:begin
          alpha[i] = -m[i]/l[i]
          beta1[i] = -n[i]/l[i]
        end
        5:begin
          alpha[i] = -l[i]/n[i]
          beta1[i] = -m[i]/n[i]
        end
      endcase
    end

; define all of the numerical constants to use in determining x and y
    r_0 = 0.577350269
    gam_s = 1.37484847732
    em = 0.004869491981
    gam = -0.13161671474
    ome = -0.159596235474
    d_0 = 0.0759196200467
    d_1 = -0.0217762490699
    c_00 = 0.141189631152
    c_10 = 0.0809701286525
    c_01 = -0.281528535557
    c_20 = -0.178251207466
    c_11 = 0.15384112876
    c_02 = 0.106959469314
    fconst = 45.0d0

    x = fconst*(alpha*gam_s+alpha^3*(1-gam_s)+alpha*beta1^2*(1-alpha^2)*$
        (gam+(em-gam)*alpha^2+(1-beta1^2)*(c_00+c_10*alpha^2+c_01*beta1^2+$
        c_20*alpha^4+c_11*alpha^2*beta1^2+c_02*beta1^4))+alpha^3*(1-alpha^2)*$
        (ome-(1-alpha^2)*(d_0+d_1*alpha^2)))
    y = fconst*(beta1*gam_s+beta1^3*(1-gam_s)+beta1*alpha^2*(1-beta1^2)*$
        (gam+(em-gam)*beta1^2+(1-alpha^2)*(c_00+c_10*beta1^2+c_01*alpha^2+$
        c_20*beta1^4+c_11*beta1^2*alpha^2+c_02*alpha^4))+beta1^3*(1-beta1^2)*$
        (ome-(1-beta1^2)*(d_0+d_1*beta1^2)))


    if noface eq 1 then begin
        xf=fconst*[0.0d0,0.0d0,2.0d0,4.0d0,6.0d0,0.0d0]
        yf=fconst*[2.0d0,0.0d0,0.0d0,0.0d0,0.0d0,-2.0d0]
        x=x+xf[face]
        y=y+yf[face]
    endif
  end
  
  'QSC':begin
; calculate direction cosines
    l = cos(theta)*sin(phi)
    m = cos(theta)*cos(phi)
    n = sin(theta)

; determine the face on which the x and y coordinates will reside by setting
; rho equal to the maximum of n,m,l,-m,-l,-n which corresponds to faces 0
; through 5 respectively
    rho = lng*0
    if size(lng,/N_dimen) EQ 0 then face = 0 else face = lonarr(n_long)

; use an array to store a remapping of the direction cosines.  This way, faces
; 0 and 5 take points on their borders with faces 1-4.  The reason for this is
; that if the max function sees identical values in an array, it takes the 
; index of the first occurrence of that value.
    remap = [0,5,2,1,4,3]

    for i = 0l, n_long-1 do begin
      dir_cos = float([n[i],-n[i],l[i],m[i],-l[i],-m[i]])
      rho[i] = max(dir_cos,temp)
      face[i] = remap[temp]
    endfor

; based on the face determined for each point, find the parameters alpha and
; beta1
    alpha = lng*0
    beta1 = alpha
    for i = 0l, n_long-1 do begin
      case face[i] of
        0:begin
          alpha[i] = l[i]/n[i]
          beta1[i] = -m[i]/n[i]
        end
        1:begin
          alpha[i] = l[i]/m[i]
          beta1[i] = n[i]/m[i]
        end
        2:begin
          alpha[i] = -m[i]/l[i]
          beta1[i] = n[i]/l[i]
        end
        3:begin
          alpha[i] = l[i]/m[i]
          beta1[i] = -n[i]/m[i]
        end
        4:begin
          alpha[i] = -m[i]/l[i]
          beta1[i] = -n[i]/l[i]
        end
        5:begin
          alpha[i] = -l[i]/n[i]
          beta1[i] = -m[i]/n[i]
        end
      endcase
    end

    x = lng*0
    y = x &  xi = y

    s = 2.d0*(((alpha gt abs(beta1)) or (beta1 ge abs(alpha))) - 0.5d0)

    case_1 = where(abs(alpha) gt abs(beta1))
    case_2 = where((abs(alpha) le abs(beta1)) and (beta1 ne 0.d0))
    case_3 = where((alpha eq 0.d0) and (beta1 eq 0.d0))
    if (case_1[0] ne -1) then xi[case_1] = beta1[case_1]/alpha[case_1]
    if (case_2[0] ne -1) then xi[case_2] = alpha[case_2]/beta1[case_2]
    if (case_3[0] ne -1) then xi[case_3] = 0.d0

    fconst=45.0d0
    u = fconst*s*sqrt((1.d0 - rho)/(1.d0 - 1.d0/sqrt(2.d0 + xi^2)))
    v = (u/1.5d1)*radeg*(atan(xi) - asin(xi/sqrt(2.d0*(1.d0 + xi^2))))
    if (case_1[0] ne -1) then begin
      x[case_1] = u[case_1]
      y[case_1] = v[case_1]
    endif
    if (case_2[0] ne -1) then begin
      x[case_2] = v[case_2]
      y[case_2] = u[case_2]
    endif
    if (case_3[0] ne -1) then begin
      x[case_3] = 0.d0
      y[case_3] = 0.d0
    endif

    if noface eq 1 then begin
        xf=fconst*[0.0d0,0.0d0,2.0d0,4.0d0,6.0d0,0.0d0]
        yf=fconst*[2.0d0,0.0d0,0.0d0,0.0d0,0.0d0,-2.0d0]
        x=(x+xf[face])
        y=(y+yf[face])
    endif
  end
  
  'TSC':begin
; calculate direction cosines
    l = cos(theta)*sin(phi)
    m = cos(theta)*cos(phi)
    n = sin(theta)

; determine the face on which the x and y coordinates will reside by setting
; rho equal to the maximum of n,m,l,-m,-l,-n which corresponds to faces 0
; through 5 respectively
    rho = lng*0
    if size(lng,/N_dimen) EQ 0 then face = 0 else face = lonarr(n_long)

; use an array to store a remapping of the direction cosines.  This way, faces
; 0 and 5 take points on their borders with faces 1-4.  The reason for this is
; that if the max function sees identical values in an array, it takes the 
; index of the first occurrence of that value.
    remap = [0,5,2,1,4,3]

    for i = 0l, n_long-1 do begin
      dir_cos = float([n[i],-n[i],l[i],m[i],-l[i],-m[i]])
      rho[i] = max(dir_cos,temp)
      face[i] = remap[temp]
    endfor

; based on the face determined for each point, find the parameters eta and xi
    eta = lng*0
    xi = eta
    for i = 0l, n_long-1 do begin
      case face[i] of
        0:begin
          eta[i] = -m[i]
          xi[i] = l[i]
        end
        1:begin
          eta[i] = n[i]
          xi[i] = l[i]
        end
        2:begin
          eta[i] = n[i]
          xi[i] = -m[i]
        end
        3:begin
          eta[i] = n[i]
          xi[i] = -l[i]
        end
        4:begin
          eta[i] = n[i]
          xi[i] = m[i]
        end
        5:begin
          eta[i] = m[i]
          xi[i] = l[i]
        end
      endcase
    endfor
    fconst = 45.0d0
    r_theta = fconst/tan(asin(rho))
    a_phi = atan(xi,-eta)
    x = r_theta*sin(a_phi)
    y = -r_theta*cos(a_phi)
    if noface eq 1 then begin
        xf=fconst*[0.0d0,0.0d0,2.0d0,4.0d0,6.0d0,0.0d0]
        yf=fconst*[2.0d0,0.0d0,0.0d0,0.0d0,0.0d0,-2.0d0]
        x=(x+xf[face])
        y=(y+yf[face])
    endif
  end
  
  'HPX':begin
    hpx_k = 3.D ; The main HEALPIX parameters (see Calabretta 2007, MNRAS)
    hpx_h = 4.D ;
    thetalim = asin((hpx_k-1)/hpx_k) 
    
    eqfaces = where( abs(theta) le thetalim, complement=polfaces)
    x=phi*0.D
    y=theta*0.D

    ; equatorial region
    if eqfaces[0] ne -1 then begin
        x[eqfaces] = phi[eqfaces]*radeg
        y[eqfaces] = 90 * hpx_k/ hpx_h * sin( theta[eqfaces])
    endif

    ;polar regions
    if polfaces[0] ne -1 then begin
       hpx_sig = sqrt ( hpx_k * (1 - abs(sin(theta[polfaces]))))
       hpx_omega = ((hpx_k mod 2 eq 1) or theta[polfaces] gt 0)*1.D
       hpx_phic = -180 + (2*floor((phi[polfaces]*radeg+180)*hpx_h/360+(1-hpx_omega)/2.)+hpx_omega)*180/hpx_h
       x[polfaces] = hpx_phic + (phi[polfaces]*radeg-hpx_phic) * hpx_sig
       y[polfaces] = 180./hpx_h *((theta[polfaces] gt 0)*2-1)* ((hpx_k+1)/2 - hpx_sig)
    endif 
  end
  else:message,strupcase(projection_type)+$
               ' is not a valid projection type.  Reset CTYPE1 and CTYPE2'
endcase

 if keyword_set(crxy) then begin
          x = x - crxy[0]
          y = y - crxy[1]
 endif

 END