wcsxy2sph.pro

basic median filter simulation

          diff2 = abs(pi2 - theta2)
          g = where((diff1 le diff2), Ng) 
          if Ng GT 0 then theta[g] = theta1[g]
          g = where( (diff2 LT diff1) , Ng) 
          if Ng GT 0 then theta[g] = theta2[g]
    endelse

  end  
  'SZP': begin

       mu = N_elements(PV2) GT 0 ? PV2[0] : 0
       phi_c = N_elements(PV2) GT 1 ? PV2[1] : 0
       theta_c = N_elements(PV2) GT 2 ? PV2[2] : 90.0

       phi_c = phi_c/radeg & theta_c = theta_c/radeg
       xp = -mu*cos(theta_c)*sin(phi_c)
       yp =  mu*cos(theta_c)*cos(phi_c)
       zp =  mu*sin(theta_c) + 1.
     
       xx = xx/radeg  &  yy = yy/radeg
       xb = (xx - xp)/zp & yb = (yy - yp)/zp
       a = xb^2 + yb^2 + 1
       b = xb*(xx - xb) + yb*(yy - yb)
       c = (xx - xb)^2 + (yy - yb)^2 - 1.
       rad = sqrt(b^2 - a*c)
       rad1 = (-b + rad)/a
       rad2 = (-b - rad)/a
       arad1 = abs(rad1)
       arad2 = abs(rad2)
       rad = rad*0.
       g = where((arad1 LE pi2) and (arad2 GT pi2), Ng )
       if Ng GT 0 then rad[g] = rad1[g]
       g = where((arad2 LE pi2) and (arad1 GT pi2), Ng )
       if Ng GT 0 then rad[g] = rad2[g]
       g = where((arad2 LE pi2) and (arad1 LE pi2), Ng )
       if Ng GT 0 then rad[g] = rad2[g] > rad1[g]
       theta = asin(rad)
       phi = atan( xx - xb*(1-sin(theta)), -(yy - yb*(1-sin(theta))) )
        end

  'TAN':begin
    sz_x = size(xx,/dimen)
    if sz_x[0] EQ 0 then theta = pi2 else $
        theta = make_array(value=pi2,dimen = sz_x)     ;Default is 90 degrees
    r = sqrt(xx^2 + yy^2)
    g = where(r GT 0, Ng)
    if Ng GT 0 then theta[g] = atan(radeg/r[g]) 
    phi = atan(xx,-yy)
  end

  'SIN':begin

    PV2_1 = N_elements(PV2) GE 1 ? PV2[0] : 0.0
    PV2_2 = N_elements(PV2) GE 2 ? PV2[1] : 0.0

    if (pv2_1 EQ 0) and (pv2_2 EQ 0) then begin
       theta = acos(sqrt(xx^2 + yy^2)/radeg)
       phi = atan(xx,-yy)
    endif else begin
       x = xx/radeg & y = yy/radeg
       a = pv2_1^2 + pv2_2^2 + 1
       b = pv2_1*(x - pv2_1) + pv2_2*(y - pv2_2)
       c = (x - pv2_1)^2 + (y - pv2_2)^2 - 1.
       rad = sqrt(b^2 - a*c)
       rad1 = (-b + rad)/a
       rad2 = (-b - rad)/a
       arad1 = abs(rad1)
       arad2 = abs(rad2)
       rad = rad*0.
       g = where((arad1 LE pi2) and (arad2 GT pi2), Ng )
       if Ng GT 0 then rad[g] = rad1[g]
       g = where((arad2 LE pi2) and (arad1 GT pi2), Ng )
       if Ng GT 0 then rad[g] = rad2[g]
       g = where((arad2 LE pi2) and (arad1 LE pi2), Ng )
       if Ng GT 0 then rad[g] = rad2[g] > rad1[g]
       theta = asin(rad)
       phi = atan( x - pv2_1*(1-sin(theta)), -(y - pv2_2*(1-sin(theta))) )
   endelse 
  end

  'STG':begin
    theta = pi2 - 2*atan(sqrt(xx^2 + yy^2)/(2.d0*radeg))
    phi = atan(xx, -yy)
  end

  'ARC':begin
    theta = pi2 - sqrt(xx^2 + yy^2)/radeg
    phi = atan(xx, -yy)
  end

  'ZPN':  begin 
   rtheta = sqrt(xx^2 + yy^2)/radeg
   phi = atan(xx, -yy)
   g = where(pv2 NE 0, Ng)
   if Ng GT 0 then np = max(g) else np =0
   pv2 = pv2[0:np]
   n = N_elements(xx)
   theta = dblarr(n)
   for i=0, n-1 do begin
      pv = pv2
      pv[0] = pv[0] - rtheta[i]
      gamma = fz_roots(pv)
; Want only the real roots
   good = where( imaginary(gamma) EQ 0, Ng) 
   if Ng EQ 0 then message,'ERROR in ZPN computation: no real roots found'
   gamma = double( gamma[good])

; If multiple real roots are found, then we seek the value closest to the 
; approximate linear solution
   
   if Ng GT 1 then begin         
        gamma1 = -pv[0]/pv[1]
        dgmin = min(abs(gamma - gamma1), dgmin_index)
        gamma = gamma[dgmin_index]   
      good = where( (gamma GE -pi2) and (gamma LE pi2), Ng)
      if Ng EQ 0 then gamma = gamma[0] else gamma = gamma[good[0]]
   endif
   theta[i] = pi2 - gamma
   endfor
  end


  'ZEA':begin
    theta = pi2 - 2.d0*asin(sqrt(xx^2 + yy^2)/(2.d0*radeg))
    phi = atan(xx,-yy)
  end

  'AIR':begin
     
    if N_elements(PV2) LT 1 then begin
      message,/informational,$
          'pv2_1 not set, using default of pv2_1 = 90 for AIR map projection'
      pv2_1 = 9.d1
    endif else pv2_1 = pv2[0]

; Numerically solve the equation for xi, by iterating the equation for xi.
; The default initial value for xi is 30 degrees, but for some values of 
; x and y, this causes an imaginary angle to result for the next iteration of
; xi.  Unfortunately, this causes the value of xi to converge to an incorrect
; value, so the initial xi is adjusted to avoid this problem.
    theta_b = pv2_1/radeg
    xi = theta_b
    zeta_b = (pi2-theta_b)/2.d0
    if (cos(zeta_b) NE 1) then $
      a = alog(cos(zeta_b))/(tan(zeta_b))^2 $
    else a = -0.5d0
    rtheta = sqrt(xx^2 + yy^2)/(2.0d*radeg) 

    repeat begin
      bad=where( abs(exp((-rtheta - a*tan(xi))*tan(xi))) gt 1)
      if (bad[0] ne -1) then xi[bad] = xi[bad]/2.d0
    endrep until (bad[0] eq -1)

    tolerance = 1.d-12
    repeat begin

      xi_old = xi
      xi = acos(exp( (-rtheta - a*tan(xi) )*tan(xi)))

    endrep until (max(abs(xi_old - xi)) lt tolerance)

;    print,rtheta,alog(cos(xi))/tan(xi) + a*tan(xi)
    theta = pi2 - 2.d0*xi
    phi = atan(xx,-yy)
  end

  'CYP':begin
    if n_elements(pv2 eq 0) then begin
      message,/informational,$
            'PV2_1 not set, using default of pv2_1 = 0 for CYP map projection'
      pv2_1 = 0.d0
    endif else pv2_1 = pv2[0]
    if N_elements(pv2) LT 2 then begin
      message,/informational,$
            'PV2_2 not set, using default of pv2_2 = 1 for CYP map projection'
      pv2_2 = 1.d0
    endif else pv2_2 = pv2[1]
    if (pv2_1 eq -pv2_2) then message,$
      'PV2_1 = -PV2_2 is not allowed for CYP map projection.'

    eta = yy/((pv2_1 + pv2_2)*radeg)
    theta = atan(eta,1) + asin(eta*pv2_1/sqrt(eta^2 + 1.d0))
    phi = xx/(pv2_2*radeg)
  end
  
  'CAR':begin
    phi = xx/radeg
    theta = yy/radeg
  end
  
  'MER':begin
    phi = xx/radeg
    theta = 2*atan(exp(yy/radeg)) - pi2
  end
  
  'CEA':begin
    if N_elements(PV2) LT 1 then message,$
      'CEA map projection requires that PV2_1 keyword be set.'
    pv2_1 = pv2[0]
    if ((pv2_1 le 0) or (pv2_1 gt 1)) then message,$
      'CEA map projection requires 0 < PV2_1 <= 1'
    phi = xx/radeg
    theta = asin(yy*pv2_1/radeg)
  end
  
  'COP':begin
    if N_elements(PV2) LT 1 then message,$
      'COP map projection requires that PV2_1 keyword be set.'
    pv2_1 =  pv2[0]
    if N_elements(PV2) LT 2 then begin
      message,/informational,$
      'PV2_2 not set, using default of PV2_2 = 0 for COP map projection'
      pv2_2=0
    endif else pv2_2 = pv2[1]
    if ((pv2_1 lt -90) or (pv2_2 gt 90) or (pv2_1 gt 90)) then message,$
 'pv2_1 and pv2_2 must satisfy -90<=PV2_1<=90, PV2_2<=90 for COP projection'
    if (pv2_1 eq -pv2_2) then message,$
 'COP projection with PV2_1=-PV2_2 is better done as a cylindrical projection'

    theta_a = pv2_1/radeg
    alpha = pv2_2/radeg
    y_0 = radeg*cos(alpha)/tan(theta_a)
    R_theta = sqrt(xx^2+(y_0-yy)^2)
    if pv2_1 LT 0 then R_theta = -R_theta
    theta = theta_a + atan(1.d0/tan(theta_a) - R_theta/$
          (radeg*cos(alpha)))
     phi = atan( xx/R_theta,(y_0-yy)/R_theta )/sin(theta_a)
  end
  
  'COD':begin
    if N_elements(pv2) LT 1 then message,$
      'COD map projection requires that PV2_1 keyword be set.'
    pv2_1 = pv2[0]
    if N_elements(pv2) LT 2 then begin
      message,/informational,$
     'PV2_2 not set, using default of PV2_2 = 0 for COD map projection'
      pv2_2 = 0
    endif else pv2_2 = pv2[1]
    if ((pv2_1 lt -90) or (pv2_2 gt 90) or (pv2_1 gt 90)) then message,$
 'pv2_1 and pv2_2 must satisfy -90<=pv2_1<=90,pv2_2<=90 for COD projection'

; use general set of equations for pv2_1 not = pv2_2
    theta_a = pv2_1/radeg

    if (pv2_2 NE 0) then begin
      alpha = pv2_2/radeg
      C = sin(theta_a)*sin(alpha)/alpha
      Y_0 = radeg*alpha/tan(alpha)/tan(theta_a)
      R_theta = sqrt(xx^2+(y_0-yy)^2)
      if pv2_1 LT 0 then R_theta = -R_theta
       theta = theta_a + alpha/(tan(alpha)*tan(theta_a))-  R_theta/radeg
; use special set of equations for pv2_1 = pv2_2
    endif else begin
      C = sin(theta_a)
      y_0 = radeg/tan(theta_a)
     R_theta = sqrt(xx^2+(y_0-yy)^2)
     if pv2_1 LT 0 then R_theta = -R_theta
      theta = theta_a + 1.0d/tan(theta_a) - R_theta/radeg
   endelse    
    phi = atan( xx/R_theta,(y_0-yy)/R_theta )/C
   end

  'COE':begin
    if N_elements(pv2) LT 1 then message,$
      'COE map projection requires that pv2_1 keyword be set.'
    pv2_1 = pv2[0]
    if N_elements(pv2) LT 2 then begin
      message,/informational,$
      'pv2_2 not set, using default of pv2_2 = 0 for COE map projection'
      pv2_2 = 0
    endif else pv2_2 = pv2[1]
    if ((pv2_1 lt -90) or (pv2_2 gt 90) or (pv2_1 gt 90)) then message,$
 'pv2_1 and pv2_2 must satisfy -90<=pv2_1<=90,pv2_2<=90 for COE projection'
    theta_a = pv2_1/radeg
    eta = pv2_2/radeg
    theta1 = (theta_a - eta)
     theta2 = (theta_a + eta)
    s_1 = sin( theta1)
    s_2 = sin( theta2)
    stheta_a = sin(theta_a)
    gamma = s_1 + s_2
    C = gamma/2
    y_0 = radeg*2.d0*sqrt(1.d0 + s_1*s_2 - gamma*stheta_a)/gamma
    R_theta = (xx^2+(y_0-yy)^2)
    if pv2_1 LT 0 then R_theta = -R_theta
    phi = 2*atan(xx/R_theta,(y_0 - yy)/R_theta)/gamma
    theta = asin((1.d0 + s_1*s_2-(xx^2+(y_0-yy)^2)*(gamma/(2.d0*radeg))^2)/gamma)
    
  end
  
  'COO':begin
    if N_elements(pv2) LT 1 then message,$
      'COO map projection requires that pv2_1 keyword be set.'
    pv2_1 = pv2[0]
    if N_elements(pv2) LT 2 then begin
      message,/informational,$
      'pv2_2 not set, using default of pv2_2 = 0 for COO map projection'
      pv2_2 = 0
    endif else  pv2_2 = pv2[1]
    if ((pv2_1 lt -90) or (pv2_2 gt 90) or (pv2_1 gt 90)) then message,$
 'pv2_1 and pv2_2 must satisfy -90<=pv2_1<=90,pv2_2<=90 for COO projection'
    theta_1 = (pv2_1 - pv2_2)/radeg
    theta_2 = (pv2_1 + pv2_2)/radeg 
    theta_a = pv2_1/radeg
 
 
; calculate value of c in simpler fashion if pv2_1 = pv2_2
    if (theta_1 eq theta_2) then c = sin(theta_1) else $
    c = alog(cos(theta_2)/cos(theta_1))/alog(tan((pi2-theta_2)/2.d0)/$
    tan((pi2-theta_1)/2.d0))

    alpha = radeg*cos(theta_1)/(c*(tan((pi2-theta_1)/2.d0))^c)
    Y_0 = alpha*(tan((pi2-theta_a)/2.d0)^c)
    R_theta = sqrt(xx^2+(y_0-yy)^2)
    if pv2_1 LT 0 then R_theta = -R_theta
     phi = atan( xx/R_theta,(y_0-yy)/R_theta )/C
    theta = pi2 - 2*atan((R_theta/alpha)^(1.d0/c))
  end
 
  'BON':begin
    if (N_elements(pv2) LT 1) then message,$
      'BON map projection requires that PV2_1 keyword be set.'
    pv2_1 = pv2[0]
    if ((pv2_1 lt -90) or (pv2_1 gt 90)) then message,$
      'pv2_1 must satisfy -90 <= pv2_1 <= 90 for BON map projection'