gdl-astrolib 2018.02.16+dfsg-1 / usr / share / gnudatalanguage / astrolib / getrot.pro

pro getrot, hdr, rot, cdelt, DEBUG = debug, SILENT = silent, ALT=alt      ;GET ROTation 
;+
; NAME:
;    GETROT
; PURPOSE:
;     Return the rotation and plate scale of an image from its FITS header
; EXPLANATION:
;     Derive the counterclockwise rotation angle, and the X and Y scale
;     factors of an image, from a FITS image header.   The input parameter 
;     may be either a FITS image header or an astrometry structure (as 
;     obtained by extast.pro)
;
; CALLING SEQUENCE:
;     GETROT, Hdr, [ Rot, CDelt, /SILENT, DEBUG =  ]   
;             or 
;     GETROT, Astr, Rot, CDelt, /SILENT, DEBUG = ]       
;
; INPUT PARAMETERS:
;     HDR - FITS Image header (string array).  Program will extract the 
;             astrometry structure
;              or
;     ASTR -  ASTROMETRY structure, of the type returned by EXTAST.
;             See the documentation for EXTAST.PRO for details.
;
; OPTIONAL OUTPUT PARAMETERS:
;       ROT - Scalar giving the counterclockwise rotation of NORTH in DEGREES 
;               from the +Y axis of the image.
;       CDELT- 2 element vector giving the scale factors in DEGREES/PIXEL in 
;               the X and Y directions.   CDELT[1] is always positive, whereas
;               CDELT[0] is negative for a normal left-handed coordinate system,
;               and positive for a right-handed system. 
;
;       If no output variables are supplied (or /DEBUG is set), then GETROT 
;       will display the rotation and plate scale at the terminal.
;
; OPTIONAL INPUT KEYWORD
; 
;       ALT - single character 'A' through 'Z' or ' ' specifying an alternate
;             astrometry system present in the FITS header.   See extast.pro
;             for more information on the ALT keyword.    Ignored if an
;             astrometry structure rather than FITS header is supplied.
;       DEBUG - if DEBUG is set, GETROT will print the rotation for both the 
;           X and Y axis when these values are unequal.  If DEBUG is set to 2, 
;           then the output parameter ROT will contain both X and Y rotations.
;
;       /SILENT - if set, then do not provide a warning about a right-handed
;           coordinate system
; PROCEDURE:
;       If the FITS header already contains CDELT (and CD or CROTA) keyword,
;       (as suggested by the Calabretta & Greisen (2002, A&A, 395, 1077) FITS 
;       standard) then this is used for the scale factor.
;       
;       If the header contains CD keywords but no CDELT keywords (as in IRAF
;       headers) then the scale factor is derived from the CD matrix.
;
;       In case of skew (different rotations of the X and Y axes), the rotations
;       are averaged together if they are less than 2 degrees.   Otherwise,
;       a warning is given and the X rotation is used.  
;
; PROCEDURES USED:
;       EXTAST, GSSS_EXTAST
; REVISION HISTORY:
;       Written W. Landsman STX January 1987 
;       Option to return both rotations added.  J. D. Offenberg, STX, Aug 1991
;       Use new astrometry structure   W. Landsman  Mar 1994
;       Recognize a GSSS header        W. Landsman  June 1994
;       Correct rotation determination with unequal CDELT values WL October 1998
;       Consistent conversion between CROTA and CD matrix  WL  October 2000
;       Correct CDELT computations for rotations near 90 deg WL November 2002
;       Preserve sign in the CDELT output  WL June 2003
;       Check if latitude/longitude reversed in CTYPE  WL  February 2004
;       Fix problem in latitude check  M.Lombardi/W.Landsman Sep 2004
;       Added ALT keyword W. Landsman May 2005
;       Account for any rotation of the native system by examining the value
;        of LONGPOLE       H. Taylor/W. Landsman
;       Account for case where X,Y rotations differ by 2*!pi WL. Aug 2011
;-
 Compile_opt IDL2
 On_error,2

 if N_params() EQ 0 then begin
        print,'Syntax: GETROT, Hdr, [ Rot, CDelt, DEBUG= , /SILENT, ALT=]'
        print,'    OR: GETROT, Astr, [ Rot, CDelt, DEBUG= , /SILENT]'
        return
 endif

 tpi = 2*!dpi
 if ~keyword_set(DEBUG) then debug = 0
 radeg = 180.0/!DPI
 sz = size(hdr,/STR)                ;Did user supply a FITS header or a CD matrix?

 if ((sz.N_dimensions eq 1) && $
     (sz.type_name EQ 'STRING')) then begin     ;FITS header?

        extast,hdr,astr,alt=alt             ;Extract astrometry from header,
        if strmid(astr.ctype[0],5,3) EQ 'GSS' then begin
                hdr1 = hdr
                gsss_stdast, hdr1
                extast, hdr1, astr
        endif
        cd = astr.cd/RADEG      ;then extract CD matrix from astrometry.
        if N_elements(cd) NE 4 then $
            message,'ERROR - Header is missing astrometry keywords CD or CDELT'

endif else $
 if ((sz.N_dimensions eq 1) and $
     (sz.type_name EQ 'STRUCT')) then begin     ;Astrometry Structure?
      astr = hdr
      cd = astr.cd/RADEG
 endif else message, $
        'ERROR - First parameter must be an image header or astrometry structure'

 if astr.cdelt[0] NE 1.0 then begin
        cdelt = astr.cdelt
        cd[0,0] *= cdelt[0] & cd[0,1] *= cdelt[0]
        cd[1,1] *= cdelt[1] & cd[1,0] *= cdelt[1]
 endif else  cd = astr.cd/RADEG                                        
 
 ctype = strmid(astr.ctype[0],0,4)
; Check if first coordinate in CTYPE is latitude
 if (ctype EQ 'DEC-') or (strmid(ctype, 1) EQ 'LAT')  then $
      cd = reverse(cd,1)
 det = cd[0,0]*cd[1,1] - cd[0,1]*cd[1,0]
 if det LT 0 then sgn = -1 else sgn = 1
 if ~keyword_set(SILENT) then if det GT 0 then $
   message,'WARNING - Astrometry is for a right-handed coordinate system',/INF
 cdelt = fltarr(2)
 if (cd[1,0] eq 0) && (cd[0,1] eq 0) then begin ;Unrotated coordinates?
   rot = 0.
   rot2 = 0.
   cdelt[0] = cd[0,0]
   cdelt[1] = cd[1,1]
 endif else begin
   rot  = atan(  sgn*cd[0,1],  sgn*cd[0,0] ) 
   rot2 = atan( -cd[1,0],  cd[1,1] )
  

 if rot NE rot2 then begin          ;Handle unequal rotations
      if keyword_set(debug) then $
        print,'X axis rotation:',rot*!RADEG, ' Y axis rotation:',rot2*!RADEG
      if debug eq 2 then rot = [rot,rot2] else begin
      
      if abs(rot - rot2)*!RADEG LT 2 then rot = (rot + rot2)/2. else $
      if abs(rot - rot2- tpi)*!radeg lt 2  then rot = (rot +rot2 -tpi)/2. else $
      if abs(rot - rot2 +tpi)*!radeg lt 2  then rot = (rot +rot2 +tpi)/2. else $
         message,/INF, $
	 'WARNING: X and Y axis rotations differ by more than 2 degrees'
      endelse
      
 endif

  cdelt[0] =   sgn*sqrt(cd[0,0]^2 + cd[0,1]^2)
  cdelt[1] =   sqrt(cd[1,1]^2 + cd[1,0]^2)
 endelse

 rot = rot*RADEG   
 if astr.longpole NE 180.0d then rot = rot + ( 180.0d - astr.longpole )
 rot = float(rot)
 cdelt = float(cdelt*RADEG)

 if N_params() EQ 1 || keyword_set(DEBUG) then begin
        if debug LT 2 then print,'Rotation (counterclockwise)',rot,' degrees'
        print,'Sampling interval X axis',cdelt[0]*3600.,' arc seconds/pixel'
        print,'                  Y axis',cdelt[1]*3600.,' arc seconds/pixel'
 endif

 return
 end