Viewing contents of file '../idllib/ssw/allpro/arcmin2hel.pro'
;+
; Project : SOHO - CDS
;
; Name : ARCMIN2HEL()
;
; Purpose : Convert arcmin from sun centre to heliographic coords.
;
; Explanation : Converts an (x,y) position given in arcmins relative to the
; solar disk centre to heliographic coordinates for the date
; supplied ( default date = today).
;
; Use : IDL> helio = arcmin2hel(xx,yy,date=date)
;
; Inputs : xx - E/W coordinate in arc minutes relative to sun center
; (West is positive); can be a vector
; yy - S/N coordinate in arc minutes relative to sun center
; (North is positive); can be a vector
;
; Opt. Inputs : None
;
; Outputs : Function returns a 2xN element vector: [lat, long] in
; degrees, where N is number of elements in XX and YY
;
; Opt. Outputs: None
;
; Keywords : date - date/time in any CDS format
;
; soho - if set uses the SOHO view point rather than
; the Earth. Note this functionality is
; duplicated by the system variable SC_VIEW,
; which in turn is set by the procedures
; USE_EARTH_VIEW or USE_SOHO_VIEW.
;
; off_limb - flag which is true if the coordinates are beyond
; the solar limb.
; error - Output keyword containing error message;
; a null string is returned if no error occurs
; radius - solar radius value [output in arcmin units]
; no_copy - input xx and yy arrays are destroyed
; p - If specified, override pb0r for p angle
; b0 - If specified, override pb0r for b angle
; r0 - If specified, override pb0r for solar apparent
; dia: r0 is the observer's distance from the
; Sun, in solar radii (one may use "zunits" to
; convert fron solar radii to, say, kilometers).
; l0 - If specified, this is the longitude (relative
; to Earth- or SOHO- central meridian) of the
; sub-observer meridian. Useful for reprojections
; and locations far from Earth.
; sphere - If set, return longitude over the whole sphere
; rather than mirror reflecting around the
; lon=90 meridian. (Useful for nonzero B angle)
;
; Calls : PB0R
; ANYTIM2UTC
;
; Restrictions: If the supplied coordinates are outside the solar disk, the
; region is projected onto the limb.
;
; Side effects: None
;
; Category : Util, coord
;
; Prev. Hist. : Original by J P Wuelser.
;
; Written : CDS version by C D Pike, RAL, 6 Sept 93
;
; Modified : Converted to use CDS time and pb0r routines, CDP, 17-May-94
; Version 3, William Thompson, GSFC, 14 November 1994
; Modified .DAY to .MJD
; Version 4, July 31, 1995, Liyun Wang, GSFC/ARC
; Vectorized the input parameters (i.e., the two
; input parameters can be vectors)
; Added ERROR keyword
; Version 5, 26-Feb-96, CDP
; Added SOHO keyword.
; Version 6, March 11, 1996, Liyun Wang, GSFC/ARC
; Modified such that point of view can be changed to
; SOHO if the env variable SC_VIEW is set to 1
; Version 7, March 13, 1996, Liyun Wang, GSFC/ARC
; Allow scalar and vector input.
; Version 8, August 2, 1996, Liyun Wang, NASA/GSFC
; Error from calling PB0R is ignored
; Version 9, June 10, 1997, Liyun Wang, NASA/GSFC
; Fixed a bug that occurs when trying to avoid off-limb points
; Version 10, Sept 10, 1997, Zarro SAC/GSFC
; Added RADIUS output keyword
; Version 11, Oct 7, 1997, Zarro SAC/GSFC
; Added some calls to TEMPORARY
; Version 12, November 4, 1997 Giulio Del Zanna, UCLAN, UK
; corrected the fact that the DATE format was changed
; inside the routine (e.g. string to a structure)
; Version 13, February 20, 1998, Zarro (SAC/GSC)
; - Added more error checking and improved memory management
; - Added /NO_COPY (use with care)
; Version 14, Dec 17, 1998, DeForest (Stanford/GSFC)
; - Added PB0R overrides on a keyword-by-keyword basis
; Version 15, 6-Jan-1999, DeForest (Stanford/GSFC)
; - Added keyword override to remove mirror-image behavior
; around the azimuthal meridian.
; Version 16, 7-Jan-99, Zarro (SMA/GSFC)
; - Made use of SOHO keyword more logical
; Version 17, 8-Apr-99, DeForest (Stanford/GSFC)
; - Added L0 keyword; renamed R -> R0 to avoid conflict with
; "radius" keyword (FDWIM).
; Version 17.5, 8-Apr-99, DeForest (Stanford/GSFC)
; - Modified interpretation of R0 keyword to agree with
; hel2arcmin: observer distance from Sun.
; Version 17.6, 9-Apr-99, Andretta (CUA/GSFC)
; - Corrected a bug in the distance-radius conversion
;
; Version : Version 17.6
;-
FUNCTION arcmin2hel, xx_in, yy_in, date=this_date, off_limb=offlimb, $
soho=soho, error=error,radius=sunr,no_copy=no_copy, $
p=p_kw,b0=b0_kw,r0=r_kw,l0=l_kw,sphere=sphere
;
; check if Soho perspective wanted.
;
use_soho=soho_view()
if exist(soho) then use_soho=soho
use_soho= 0 > use_soho < 1
error = ''
dummy=[-99., -99.]
IF N_ELEMENTS(xx_in) NE N_ELEMENTS(yy_in) THEN BEGIN
error = 'Two input parameter not compatible'
MESSAGE, error, /cont
RETURN, dummy
ENDIF
;-- use temporary to save memory I/O
if keyword_set(no_copy) then begin
xx = temporary(xx_in)
yy = temporary(yy_in)
endif else begin
xx=xx_in
yy=yy_in
endelse
xs = SIZE(xx)
ys = SIZE(yy)
IF xs(0) GT 1 THEN BEGIN
IF xs(1) NE 1 THEN BEGIN
error = 'Parameter must be a column or row vector.'
MESSAGE, error, /cont
RETURN,dummy
ENDIF ELSE xx = TRANSPOSE(temporary(xx))
ENDIF
IF ys(0) GT 1 THEN BEGIN
IF ys(1) NE 1 THEN BEGIN
error = 'Parameter must be a column or row vector.'
MESSAGE, error, /cont
RETURN, dummy
ENDIF ELSE yy = TRANSPOSE(temporary(yy))
ENDIF
offlimb = INTARR(N_ELEMENTS(xx))
;----------------------------------------------------------------------
; Deal with P, B0, and R. If all three are specified, don't bother
; finding out where the Earth was...
;----------------------------------------------------------------------
if not (isvalid(p_kw) and isvalid(b0_kw) and isvalid(r_kw)) then begin
;---------------------------------------------------------------------
; set up date required in format required by PB0R
;---------------------------------------------------------------------
IF N_ELEMENTS(this_date) EQ 0 THEN get_utc, dat ELSE $
dat = anytim2utc(this_date)
IF dat.mjd EQ 0 THEN BEGIN
error = 'Wrong time format specification.'
MESSAGE, error, /cont
RETURN, dummy
ENDIF
;--------------------------------------------------------------------
; get B0 and solar radius
;--------------------------------------------------------------------
angles = float(pb0r(dat, soho=use_soho, error=error))
;--------------------------------------------------------------------
; normalize, check if off limb, project back to limb
;--------------------------------------------------------------------
if isvalid(p_kw) then angles(0) = p_kw
if isvalid(b0_kw) then angles(1) = b0_kw
if isvalid(r_kw) then angles(2) = !radeg*60 * atan(1.0/r_kw); radians-to-arcmin conversion
end else angles = [p_kw,b0_kw,!radeg*60 * atan(1.0/r_kw)]
sunr = angles(2)
b0 = angles(1)/!radeg
xx = temporary(xx/sunr) & yy = temporary(yy/sunr)
inr = SQRT(xx*xx+yy*yy)
ii = WHERE(inr GT 1.0)
IF ii(0) GE 0 THEN BEGIN
xx(ii) = xx(ii)/inr(ii)
yy(ii) = yy(ii)/inr(ii)
inr(ii) = 1.0
offlimb(ii) = 1
ENDIF
;---------------------------------------------------------------------------
; assume inxy on sphere, calc. x,y,z. Here: x towards obs., y = W, z = N
;---------------------------------------------------------------------------
xyz = TRANSPOSE([[SQRT(1.0-inr*inr)], [xx], [yy] ])
;---------------------------------------------------------------------------
; rotate around y axis to correct for B0 angle (B0: hel. lat. of diskcenter)
;---------------------------------------------------------------------------
rotmx = [[COS(b0), 0.0, SIN(b0)], [0.0, 1.0, 0.0], [-SIN(b0), 0.0, COS(b0)]]
xyz = rotmx # temporary(xyz)
xyz = TRANSPOSE(temporary(xyz))
;---------------------------------------------------------------------------
; calc. latitude and longitude.
;---------------------------------------------------------------------------
latitude = ASIN(xyz(*,2))
latitude = temporary(latitude) < (89.99/!radeg) ; force lat. smaller 89.99 deg.
latitude = temporary(latitude) > (-89.99/!radeg) ; force lat. bigger -89.99 deg.
longitude = ATAN(xyz(*,1), xyz(*,0)) ; longitude
sphere=1
if not keyword_set(sphere) then begin
;---------------------------------------------------------------------------
; longitude may be larger than 90 degrees due to nonzero B0: get proper value
;---------------------------------------------------------------------------
ii = WHERE(xyz(*, 0) LT 0.0)
IF ii(0) GE 0 THEN BEGIN
tmp = xyz(ii, *)
tmp_l = longitude(ii)
jj = WHERE(tmp(*, 1) GE 0.0)
IF jj(0) GE 0 THEN BEGIN
tmp_l(jj) = !pi-tmp_l(jj)
longitude(ii) = tmp_l
ENDIF
tmp_l = longitude(ii)
jj = WHERE(tmp(*, 1) LT 0.0)
IF jj(0) GE 0 THEN BEGIN
tmp_l(jj) = -!pi-tmp_l(jj)
longitude(ii) = tmp_l
ENDIF
ENDIF
end
;---------------------------------------------------------------------------
; convert to degrees
;---------------------------------------------------------------------------
delvarx,xyz,inr,xx,yy
if not isvalid(l_kw) then $
RETURN, !radeg*TRANSPOSE([[latitude], [longitude]]) $
else $
RETURN, !radeg*TRANSPOSE([[latitude], [longitude+l_kw]])
END