C C This sample program uses the FITSIO subroutines to add World Coordinate C System (WCS) information to a full-disk solar image. C C Author: C C Stephen Walton C Department of Physics and Astronomy C California State University Northridge C 18111 Nordhoff St. C Northridge, CA 91330-8268 C stephen.walton@csun.edu C C Version 1.0. C C This program has been divided into two parts. A main program simply C opens a FITS file and calls a subroutine ADDWCS. This subroutine C actually adds the WCS information to the file. Exactly what it adds C and why is documented in the large comment block at the start of that C subroutine. C C Syntactically (and I hope semantically), this program is strictly standard C Fortran-77 except for the two IMPLICIT NONE statements. C PROGRAM WCSSMP C C The necessary variables. The IMPLICIT NONE statement below is not C standard Fortran-77, so I've commented it out (though I left it in when C I tested this code). C IMPLICIT NONE INTEGER STATUS, UNIT, RWMODE, BSIZE CHARACTER*80 FILNAM, ERRMSG C C Open the FITS file. This portable program prompts the user for a file C name, but on many systems the below could be replaced by C CALL GETARG(1, FILNAM) C 1 WRITE (*, *) 'Enter file name' READ (*, '(A80)') FILNAM IF (FILNAM .EQ. ' ') THEN WRITE (*, *) 'Try again.' GOTO 1 ENDIF C C FITSIO calls abort if the STATUS argument is nonzero, so we must C initialize it to zero. We set the RWMODE argument to 1 to C indicate that the file is to be opened with write permission. C Then, we find an unused I/O unit and open the file. C STATUS=0 RWMODE=1 CALL FTGIOU(UNIT, STATUS) CALL FTOPEN(UNIT, FILNAM, RWMODE, BSIZE, STATUS) IF (STATUS .NE. 0) THEN WRITE (*, *) 'Fatal error: Could not open file ', FILNAM STOP ENDIF CALL ADDWCS(UNIT, STATUS) C C If an error occurs in ADDWCS, it will return the FITSIO error value C in the STATUS variable. We get the corresponding message, then C clear STATUS so that we can close the file. C IF (STATUS .NE. 0) THEN CALL FTGMSG(ERRMSG) WRITE (*, *) 'FATAL ERROR NUMBER', STATUS WRITE (*, *) ERRMSG STATUS = 0 ENDIF C C Call the subroutine which actually adds the WCS values to the file. C C C Almost done now. Close the image and free the I/O unit. C CALL FTCLOS(UNIT, STATUS) CALL FTFIOU(UNIT, STATUS) STOP END SUBROUTINE ADDWCS(UNIT, STATUS) C C This subroutine adds World Coordinate System (WCS) information to C the FITS file open on FITSIO unit UNIT, assuming that image is a C Kitt Peak Vacuum Tower telescope image. C C C First, I'll start with a description of exactly which assumptions this code C makes. The image to which a World Coordinate System is to be added is C assumed to be one from the Kitt Peak Vacuum Tower Telescope. These FITS C files are images such that the first row stored is the southernmost one, and C the last pixel in each row is the westernmost one. In other words, if the C image is displayed with the first datum in the image in the bottom left C corner of the display, the displayed image will have solar north up and C solar west to the right. C C The geometry of the image is described by five keywords, produced by the C FNDLMB program (hence their names), which give the size of an ellipse fit C to the limb of the image: C C FNDLMBMI--the length of the minor axis (pixels) C FNDLMBMA--the length of the major axis (pixels) C FNDLMBXC--the X (column) coordinate of the image center (pixels) C FNDLMBYC--the Y (row) coordinate of the image center (pixels) C FNDLMBAN--the angle between the image minor axis and the row C direction (degrees) C C With this information, we can translate an image pixel point (ix, iy) into a C point in normalized solar coordinates (x, y) defined such that x**2 + y**2=1 C for points on the limb of the Sun, (x=1, y=0) is the north limb point, and C (x=0, y=1) is the west limb point. I find it helpful to think of this C transformation in three steps: First, rotate the coordinate axes by an angle C FNDLMBAN to translate (ix, iy) into (ix', iy'). These two points now satisfy C the equation (ix'/FNDLMBMI)**2 + (iy'/FNDLMBMA)**2=1. The second step is to C produce (x', y') from (ix', iy') by dividing ix' by FNDLMBMI and iy' by C FNDLMBMA. Then (x'**2 + y'**2)=1. Finally, the coordinate axes are rotated C by an angle -FNDLMBAN for the final translation from (x', y') to (x, y). C As a series of matrices: C C ( x ) ( cos(AN) -sin(AN) )( 1/MI 0 )( cos(AN) sin(AN) )( ix ) C ( ) = ( )( )( )( ) C ( y ) ( sin(AN) cos(AN) )( 0 1/MA )(-sin(AN) cos(AN) )( iy ) C C where, for the sake of brevity, the leading FNDLMB has been removed from the C keyword names. Notice that here (ix=i-FNDLMBXC, iy=j-FNDLMBYC) are the C offsets from image center. Carrying out the above matrix multiplication C gives the matrix elements which are written to the image header in the code C below. C IMPLICIT NONE INTEGER UNIT, STATUS C C Declare local variables C DOUBLE PRECISION XC, YC, MA, MI, AN, PI, PC(2, 2), C, S CHARACTER*80 COMM C C Read the header keyword values. As noted above, the STATUS value, C if set to a non-zero error by any of the calls below, will not be C changed, and all subsequent calls will return without doing anything. C Thus we only check the STATUS after the group is completed. C CALL FTGKYD(UNIT, 'FNDLMBXC', XC, COMM, STATUS) CALL FTGKYD(UNIT, 'FNDLMBYC', YC, COMM, STATUS) CALL FTGKYD(UNIT, 'FNDLMBMA', MA, COMM, STATUS) CALL FTGKYD(UNIT, 'FNDLMBMI', MI, COMM, STATUS) IF (STATUS .NE. 0) RETURN C C For maximum flexibility, we allow the ANGLE argument to be missing, C in which case it is assumed to be zero. C CALL FTGKYD(UNIT, 'FNDLMBAN', AN, COMM, STATUS) IF (STATUS .NE. 0) THEN AN = 0.0D0 STATUS = 0 ENDIF C C Calculate the four matrix elements above. C PI = 4.0D0*ATAN(1.0D0) C = COS(PI*AN/180.0D0) S = SIN(PI*AN/180.0D0) PC(1, 1) = C**2/MI + S**2/MA PC(1, 2) = C*S*(1.0D0/MI - 1.0D0/MA) PC(2, 1) = PC(1, 2) PC(2, 2) = C**2/MA + S**2/MI C C PORTABILITY NOTE: We write the PC matrix as the older CD matrix standard. C This is to maintain compatibility with IRAF and other older software. C Once the current draft WCS standard is made final and your software is C updated to read the new format, you can change the CD matrix below to C the PC matrix by changing 'CD1_1' to 'PC001001' and so on. C CALL FTPKYD(UNIT, 'CD1_1', PC(1, 1), 10, ' ', STATUS) CALL FTPKYD(UNIT, 'CD1_2', PC(1, 2), 10, ' ', STATUS) CALL FTPKYD(UNIT, 'CD2_1', PC(2, 1), 10, ' ', STATUS) CALL FTPKYD(UNIT, 'CD2_2', PC(2, 2), 10, ' ', STATUS) IF (STATUS .NE. 0) RETURN C C Now we copy the center coordinates into the keywords used for WCS. C CALL FTGKYD(UNIT, 'FNDLMBXC', XC, COMM, STATUS) CALL FTPKYD(UNIT, 'CRPIX1', XC, 10, COMM, STATUS) CALL FTGKYD(UNIT, 'FNDLMBYC', YC, COMM, STATUS) CALL FTPKYD(UNIT, 'CRPIX2', YC, 10, COMM, STATUS) C C This ends the addition of the linear part of the WCS. Additional C values can be written here to define a non-linear WCS which would C map the image values to solar longitude and latitude. Since this C depends on a good deal of software support which does not yet C exist, I will only outline what should be done without doing it. C C A proposal is before the FITS Working Group to define PLON and PLAT C as planetary (and solar, but SLON/SLAT were taken) longitude and C latitude, respectively. For a spherical object viewed from a C finite distance, the map projection is known as azimuthal (AZP). C To add this information to the image, one would create CTYPE1 and C CTYPE2 keywords with values 'PLON-AZP' and 'PLAT-AZP', respectively. C CRVAL1 and CRVAL2 would be, respectively, the longitude and latitude C of the pixel whose location is in CRPIX1 and CRPIX2. For a typical C full-disk solar image, these would be the sub-earth (or sub-spacecraft) C latitude and longitude. The special keyword PROJP1 would give the C distance of the observer from the center of the Sun in units of the C solar radius; it would be negative to indicate the observer is looking C at a sphere from the outside, as opposed to the "usual" case of the C celestial sphere observed from the inside. C RETURN END