PRO wcsxy2sph, x, y, longitude, latitude, map_type, ctype = ctype, $
face=face, pv2 = pv2, pv1 = pv1, $
crval=crval, crxy = crxy, longpole = longpole, Latpole = latpole
compile_opt idl2, hidden
; Define angle constants
pi = !DPI
radeg = 180d0/!dpi
pi2 = pi/2.0d
;
; Please keep the following list up to date in extast as well as here
; and in wcssph2xy:
;
map_types=['DEF','AZP','TAN','SIN','STG','ARC','ZPN','ZEA','AIR','CYP',$
'CAR','MER','CEA','COP','COD','COE','COO','BON','PCO','SFL',$
'PAR','AIT','MOL','CSC','QSC','TSC','SZP','HPX','HCT','XPH']
; check to see that enough parameters (at least 4) were sent
if ( N_params() lt 4 ) then begin
print,'Syntax - WCSXY2SPH, x, y, longitude, latitude,[ map_type, '
print,' CTYPE= , FACE=, PV2= , CRVAL= , CRXY= , '
print,' LATPOLE = , LONGPOLE = ]'
return
endif else if (n_params() eq 5) then begin
if keyword_set(ctype) then message,$
'Use either the MAP_TYPE positional parameter or set the projection type' + $
'with CTYPE, but not both.'
; set projection_type string using map_type parameter (a number)
ntypes = n_elements(map_types)
if (N_ELEMENTS(map_type) eq 1 && map_type ge 0 && $
map_type lt ntypes) then begin
projection_type = map_types[map_type]
endif else message,'MAP_TYPE must be a scalar >= 0 and < '+$
strtrim(string(ntypes),2)+'; it was set to '+$
strtrim(string(map_type),2)
endif else if (n_params() eq 4) then wcs_check_ctype, ctype, projection_type
; checks CTYPE format and extract projection type
; GENERAL ERROR CHECKING
; find the number of elements in each of the data arrays
n_x = n_elements(x)
n_y = n_elements(y)
sz_x = size(x)
sz_y = size(y)
; check to see that the data arrays have the same size
if (n_x ne n_y) then $
message,'The arrays X and Y must have the same number of elements.'
; this sets the default map projection type for the cases when map_type or
; projection_type is set to 'DEF' or if projection_type is not set at this
; point. As suggested in 'Representations of Celestial Coordinates in FITS'
; the default type is set to CAR (Plate Carree) the simplest of all projections.
if (n_elements(projection_type) eq 0) || (projection_type eq 'DEF') then $
projection_type='CAR'
; Check to make sure all the correct parameters and keywords are set for
; spherical projections.
if ((N_elements(ctype) EQ 3) || keyword_set(face) || $
(projection_type eq 'CSC') || $
(projection_type eq 'QSC') || (projection_type eq 'TSC')) then begin
noface = ~keyword_set(face)
endif
; check to see if the x and y offsets are set properly. If not, break out
; of program. If so, apply offsets. If the x and y offsets are not set,
; then assume they are zero.
if ( (keyword_set(crxy)) && N_elements(crxy) NE 2) then $
message,'Offset keyword CRXY must contain 2 elements'
if keyword_set(crxy) && ~array_equal(crxy,[0d0,0d0]) then begin
xx = double(x + crxy[0] )
yy = double(y + crxy[1] )
endif else begin
xx = double(x)
yy = double(y)
endelse
if ( N_elements(crval) eq 1 ) then $
message,'CRVAL keyword must be a 2 element vector'
; BRANCH BY MAP PROJECTION TYPE
case strupcase(projection_type) of
'AZP':begin
PV2_1 = N_elements(PV2) GE 1 ? PV2[0] : 0.0 ; PV2_1 =mu (spherical radii)
PV2_2 = N_elements(PV2) GE 2 ? PV2[1] : 0.0 ; PV2_2 = gamma (degrees)
if (pv2_1 lt 0) then message,$
'AZP map projection requires the keyword pv2_1 >= 0'
gamma = pv2_2/radeg
mu = pv2_1
r = sqrt(xx^2 + yy^2*cos(gamma)^2)
rho = r/(radeg*(mu+1) + yy*sin(gamma) )
omega = asin( rho*mu/ sqrt( rho^2 + 1.d0) )
xsi = atan(1.d0, rho)
phi = atan(xx, -yy*cos(gamma) )
theta1 = xsi - omega
theta2 = xsi + omega + !dpi
theta = theta1*0.0
if abs(mu) LT 1 then begin
g = where(abs(theta1) LT pi2, Ng)
if Ng GT 0 then theta[g] = theta1[g]
g = where(abs(theta2) LT pi2, Ng)
if Ng GT 0 then theta[g] = theta2[g]
endif else begin
diff1 = abs(pi2 - theta1)
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) && (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
if size(yy,/N_dimen) EQ 0 then theta = theta[0] ;Make scalar again
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) || (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) || (pv2_2 gt 90) || (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) || (pv2_2 gt 90) || (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) || (pv2_2 gt 90) || (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) || (pv2_2 gt 90) || (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) || (pv2_1 gt 90)) then message,$
'pv2_1 must satisfy -90 <= pv2_1 <= 90 for BON map projection'
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 initial 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) || $
(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
'GLS':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]
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 ; HEALPix projection
; See Calabretta & Roukema 2007, MNRAS, 381, 865
pv2_1 = N_ELEMENTS(pv2) GE 1 ? pv2[0] : 4.d
pv2_2 = N_ELEMENTS(pv2) GE 2 ? pv2[1] : 3.d
hpx_k = pv2_2 ; The main generalised HEALPIX parameters
hpx_h = pv2_1 ;
ik = ROUND(hpx_k)
ih = ROUND(hpx_h)
IF (ik le 0 || ih le 0) THEN MESSAGE, $
'Illegal PV2 array:' + STRCOMPRESS(STRJOIN(pv2,/single))+ $
'; should be positive integers'
phi = xx ; Create theta & phi arrays in same shape as xx & yy.
theta = yy ;
; Is pixel is in an unoccupied facet ?
invfd = ih / 360d0 ; inverse facet diagonal
wdiag = ik / 2 ; Semi-width of occupied diagonal in facets
IF ik THEN BEGIN ; Odd K (including standard HEALPix):
xoff = ih + wdiag + 1 ? 0.5d0 : 0d0
ix = ROUND( (xx + yy)*invfd + xoff ) ; Facet indices
diag = ROUND( (yy - xx)*invfd - xoff ) + wdiag + TEMPORARY(ix)
ENDIF ELSE BEGIN ; Even K
; Row offset for (0,0) depends on
; whether h is odd or even:
yoff = ih + (ik-1)/2 ? -0.25d0 : 0.25d0
ioff = ih ? 2*((ik-2)/4) + 1 : 2*(ik/4)
ix = ROUND( (xx + yy)*invfd + yoff )
diag = ROUND( (yy - xx)*invfd + yoff ) + ioff + TEMPORARY(ix)
ENDELSE
hpx_good = WHERE(diag LT ik AND diag GE 0, COMPLEMENT = hpx_bad)
if hpx_bad[0] ne -1 then begin ; Set coords of off-sky pixels to NaN:
phi[hpx_bad] = !VALUES.D_NAN
theta[hpx_bad] = !VALUES.D_NAN
endif
ylim = 90*(ik-1)/hpx_h
eqfaces = -1
polfaces = -1
IF hpx_good[0] NE -1 THEN BEGIN
equas = where(abs(yy[hpx_good]) le ylim, complement=poles)
IF equas[0] NE -1 THEN eqfaces = hpx_good[TEMPORARY(equas)]
IF poles[0] NE -1 THEN polfaces = hpx_good[TEMPORARY(poles)]
ENDIF
hpx_good = 0
; equatorial region
if eqfaces[0] ne -1 then begin
phi[eqfaces]=xx[eqfaces]/radeg
theta[eqfaces]=asin(yy[eqfaces]*(hpx_h/(ik*90)))
; Allow wrapped values of x, so following commented out:
; hpx_bad = where(xx[eqfaces] lt -180.D or xx[eqfaces] gt 180.D)
; if hpx_bad[0] ne -1 then begin
; phi[eqfaces[hpx_bad]]=!VALUES.D_NAN
; theta[eqfaces[hpx_bad]]=!VALUES.D_NAN
; endif
endif
; polar regions
if polfaces[0] ne -1 then begin
hpx_sig = (ik+1)/2.D - abs(yy[polfaces])*(hpx_h/180d)
hpx_omega = FIX((hpx_k mod 2 eq 1) or yy[polfaces] gt 0)
hpx_xc = -180 + (2 * floor( (xx[polfaces]+180d0)*hpx_h/360d0 + $
(1-hpx_omega)/2d0 ) $
+ hpx_omega) * 180d0/hpx_h
poles = where(hpx_sig EQ 0.D) ; Avoid divide by zero at poles
IF poles[0] NE -1 THEN hpx_sig[poles] = 1.D
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)
IF poles[0] NE -1 THEN BEGIN
phi[polfaces[poles]] = hpx_xc[poles] / radeg
theta[polfaces[poles]] = $
((yy[polfaces[poles]] gt 0)*2-1) * !dpi/2d0
ENDIF
endif
end
'HCT':begin
phi = xx/radeg
theta = dblarr(n_elements(yy))
ylim = 90*(3-1)/4
w_np = where(yy ge ylim, n_np)
w_eq = where((yy lt ylim) and (yy gt -ylim), n_eq)
w_sp = where(yy le -ylim, n_sp)
if n_np gt 0 then theta[w_np] = asin(1-(2-yy[w_np]/ylim)^2/3.d)
if n_eq gt 0 then theta[w_eq] = asin((yy[w_eq]/ylim)*2./3.d)
if n_sp gt 0 then theta[w_sp] = -asin(1-(2+yy[w_sp]/ylim)^2/3.d)
end
'XPH':begin ; Butterfly re-arrangement of HPX, see Calabretta & Lowe (2013)
;
; xx, yy reference point is the pole which is in the centre of the
; grid. Diagonal length in IWC is 360deg/sqrt(2)
;
; identify on-sky & off-sky pixels:
scale = 1d0/sqrt(2d0)
halfwidth = 180d0*scale
quarter = 90d0*scale
dg1 = (xx + yy)
dg2 = (xx - yy)
good = WHERE( xx ge -halfwidth AND xx le halfwidth AND $
yy ge -halfwidth AND yy le halfwidth AND $
((dg1 ge -quarter AND dg1 le quarter) OR $
(dg2 ge -quarter AND dg2 le quarter)), ngood, $
COMPLEMENT = bad, NCOMPLEMENT = nbad)
scalar = sz_x[0] EQ 0
dims = scalar ? 1 : sz_x[1:sz_x[0]]
phi = DBLARR(dims,/NOZERO) & theta = DBLARR(dims,/NOZERO)
IF nbad GT 0 THEN BEGIN
phi[bad] = !values.D_NAN
theta[bad] = !values.D_NAN
bad = 0
ENDIF
IF ngood GT 0 THEN BEGIN
dg1 = dg1[good]*scale
dg2 = dg2[good]*scale
xi = DBLARR(ngood,/NOZERO)
eta = DBLARR(ngood,/NOZERO)
phig = DBLARR(ngood)
; transform xy coords from XPH to HPX
xq = xx[good] ge 0d0 ; true on LHS of projection
yq = yy[good] ge 0d0 ; true in top half of projection
quad = WHERE(~xq AND yq) ; upper right
IF quad[0] ne -1 THEN BEGIN
xi[quad] = -dg1[quad]
eta[quad] = dg2[quad]
phig[quad] = -180d0
ENDIF
quad = WHERE(~xq AND ~yq) ; lower right
IF quad[0] ne -1 THEN BEGIN
xi[quad] = dg2[quad]
eta[quad] = dg1[quad]
phig[quad] = -90d0
ENDIF
quad = WHERE(xq AND ~yq) ; lower left
IF quad[0] ne -1 THEN BEGIN
xi[quad] = dg1[quad]
eta[quad] = -dg2[quad]
phig[quad] = 0d0
ENDIF
quad = WHERE(xq AND yq) ; upper left
IF quad[0] ne -1 THEN BEGIN
xi[quad] = -dg2[quad]
eta[quad] = -dg1[quad]
phig[quad] = 90d0
ENDIF
xq = 0 & yq = 0 & dg1 = 0 & dg2 = 0 & quad = 0
eta += 90d0
thetag = DBLARR(ngood,/NOZERO)
poles = WHERE(ABS(eta) gt 45d0, COMPLEMENT=equas)
IF equas[0] NE -1 THEN BEGIN
phig[equas] += xi[equas]
thetag[equas] = ASIN(eta[equas]*(2d0/135d0))
ENDIF
equas = 0
IF poles[0] NE -1 THEN BEGIN
hpx_sigma = (90d0 - ABS(eta[poles]))/45d0
test = WHERE(hpx_sigma lt 1d-4,ntest)
phig[poles] += xi[poles]/hpx_sigma
sgn = 2*(eta[poles] gt 0d0) - 1
thetag[poles] = TEMPORARY(sgn)*ASIN(1d0-hpx_sigma^2/3d0)
IF ntest GT 0 THEN $
thetag[poles[test]] = pi2 - SQRT(2d0/3d0)*hpx_sigma[test]
ENDIF
phig += 45d0
poles = 0 & hpx_sigma = 0
phi[good] = TEMPORARY(phig)/radeg
theta[good] = TEMPORARY(thetag)
good = 0
ENDIF
IF scalar THEN BEGIN
phi = phi[0]
theta = theta[0]
ENDIF
end
else:message,strupcase(projection_type) + $
' is not a valid projection type. Reset CTYPE'
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(map_type) EQ 0 then $
map_type = where(projection_type EQ map_types)
map_type = map_type[0]
conic = (map_type GE 13) && (map_type LE 16)
zenithal = ((map_type GE 1) && (map_type LE 8)) || $
(map_type EQ 26) || (map_type EQ 29)
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, pv1=pv1, /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
good = WHERE(FINITE(longitude), ngood)
IF ngood GT 0 THEN BEGIN
lgood = longitude[good]
temp = where(lgood lt 0.d0, Nneg)
if (Nneg GT 0) then lgood[temp] = lgood[temp] + 3.6d2
temp = where(lgood ge 3.6d2, Nneg)
if (Nneg GT 0) then lgood[temp] = lgood[temp] - 3.6d2
longitude[good] = lgood
ENDIF
return
end