}A CGIS-FOR.BCK CGIS-FOR.BCK,BACKUP/BLOCK=2048 [...] [-]CGIS-FOR.BCK/SAVE WEILAND  hn6WV6.2 _BORNEO::  _BORNEO$DKB400: V6.2  w*[COBE_FOR]AACGIS_FOR.TXT;1+,7 . / h 4 - 0123KPWO 562+72+89G hHJCOBE Guest Investigator Software - FORTRAN Library (Standalone) DATA I/O ________________ This library contains eight examples of FORTRAN routines for reading FITS files. These programs use the FITSIO subroutines developed by the High Energy Astrophysics Science Archive Research Center (HEASARC) at NASA Goddard Space Flight Center. The purpose of these routines are to act as a guide to a more complete FORTRAN Data I/O routine. There is one routine for the DMR instrument (DMR_READ.FOR), two routines for FIRAS (FIRAS_READ.FOR and FITS_READ_COVAR.FOR) and four routines for the DIRBE instrument, namely DIRBE_GP_READ.FOR, DIRBE_WK_READ.FOR, DIRBE_E90_READ.FOR, and SKYINFO.FOR. The differing DIRBE routines reflect the wide varying nature of the particular data products, i.e. the Galactic Plane Maps, Weekly Averaged Maps, Solar Elongation 90 maps, and the SKYMAP_INFO primary data array. The eighth routine, RASTER_READ.FOR, will read simple FITS images as opposed to the binary table format used in the COB> CGIS-FOR.BCK7  [COBE_FOR]AACGIS_FOR.TXT;1 E project data sets and initial data products. The FITSIO library is available by anonymous ftp to legacy.gsfc.nasa.gov in the /fits_info/software directory. See documentation accompanying the FITSIO library for compiliation instructions for your particular operating system. Two examples (UNIX and VMS) follow: UNIX> f77 -o dirbe_wk_read dirbe_wk_read.for -lm -L/usr/local/fitsio -lfitsio VMS> fortran/ext dirbe_wk_read.for VMS> link dirbe_wk_read,[fitsio_path]fitsio/lib An additional routine (READ_FSS.FOR) is provided as an example for reading the FIRAS native format (VAX binary) FSS_SSSKY. COORDINATE CONVERSIONS ______________________ There are two main coordinate conversion routines, PIX2LL.FOR and LL2PIX.FOR. PIX2LL returns longitude and latitude pointing to center of pixel given pixel number and resolution of the cube. Additionally the spherical quadralateralized cube face ,x, y can be returned. LL2PIX returns the pixel number, for a specified skycube resolution, given a latitude and longitude. The coordinate system is an input for both routines. Note, all coordintes are at epoch=J2000. Both routines have a number of subroutines. To aid in the compilation example VMS command files (PIX.COM and LL.COM) are given and example UNIX make files (MAKE.PIX and MAKE.LL) are supplied. In both cases the executable will be called pix2ll.exe and ll2pix.exe, respectively. As a further tool the stand alone program PIXEL_XY is supplied. This will compute the COBE quadralater:S CGIS-FOR.BCK7  [COBE_FOR]AACGIS_FOR.TXT;1 alized cube face, x, and y positions for all the DIRBE pixels. This is equivalent to reading the primary array in the DIRBE_SKYMAP_INFO.FITS file. While these routines are fully executable FORTRAN programs, the purpose of these routines are to act as a starting point for the researchers particular needs. Your comments and suggestions are welcome. If you report installation difficulties or software bugs, other users will benefit from your experience. For assistance, please contact David Leisawitz at leisawitz@stars.gsfc.nasa.gov. *[COBE_FOR]AXISXY.FOR;1+,8 . / h 4 R- 0123KPWO 568878889G hHJ{ CGIS-FOR.BCK8  [COBE_FOR]AXISXY.FOR;1 at SUBROUTINE AXISXY(C,NFACE,X,Y) C C CONVERTS UNIT VECTOR C INTO NFACE NUMBER (0-5) AND X,Y C IN RANGE 0-1 C REAL*4 C(3) AC3=ABS(C(3)) AC2=ABS(C(2)) AC1=ABS(C(1)) IF(AC3.GT.AC2) THEN IF(AC3.GT.AC1) THEN IF(C(3).GT.0.) GO TO 100 GO TO 150 ELSE IF(C(1).GT.0.) GO TO 110 GO TO 130 ENDIF ELSE IF(AC2.GT.AC1) THEN IF(C(2).GT.0.) GO TO 120 GO TO 140 ELSE IF(C(1).GT.0.) GO TO 110 GO TO 130 ENDIF ENDIF 100 NFACE=0 ETA=-C(1)/C(3) XI=C(2)/C(3) GO TO 200 110 NFACE=1 XI=C(2)/C(1) ETA=C(3)/C(1) GO TO 200 120 NFACE=2 ETA=C(3)/C(2) XI=-C(1)/C(2) GO TO 200 130 NFACE=3 ETA=-C(3)/C(1) XI=C(2)/C(1) GO TO 200 140 NFACE=4 ETA=-C(3)/C(2) XI=-C(1)/C(2) GO TO 200 150 NFACE=5 ETA=-C(1)/C(3) XI=-C(2)/C(3) GO TO 200 C 200 CALL INCUBE(XI,ETA,X,Y) X=(X+1.)/2. Y=(Y+1.)/2. RETURN END *[COBE_FOR]BIT_SET.FOR;1+,9 . / h 4 - 0123KPWO 56J287J2889G hHJbTU CGIS-FOR.BCK9  [COBE_FOR]BIT_SET.FOR;1  Subroutine BIT_SET(IX,IY,LENGTH) C C Routine to set up the bit tables for use in the pixelization subroutines C (extracted and generalized from existing routines). C C Arguments: C Integer*2 LENGTH ! Number of elements in IX, IY Integer*4 IX(length) ! X bits Integer*4 IY(length) ! Y bits C C Variables: C Integer*4 i,j,k,ip,id ! Loop variables C Do 30 I = 1,length j = i - 1 k = 0 ip = 1 10 if (j .eq. 0) then ix(i) = k iy(i) = 2*k else id = mod(j,2) j = j/2 k = ip * id + k ip = ip * 4 go to 10 endif 30 continue C return end *[COBE_FOR]CONV_E2G.FOR;1+,: . / h 4 - 0123KPWO 56%*57%*589G hHJsV CGIS-FOR.BCK:  [COBE_FOR]CONV_E2G.FOR;1  Subroutine conv_E2G(ivector,ovector) C C Routine to convert ecliptic (celestial) coordinates at the given C epoch to galactic coordinates. C C -------------- C C Galactic Coordinate System Definition (from Zombeck, "Handbook of C Space Astronomy and Astrophysics", page 71): C C North Pole: 12:49 hours right ascension (1950.0) C +27.4 degrees declination (1950.0) C C Reference Point: 17:42:26.603 hours right ascension (1950.0) C -28 55' 00.45 degrees declination (1950.0) C C -------------- C C Implicit NONE C C Arguments: C Real*4 IVECTOR(3) ! Input coordinate vector Real*4 OVECTOR(3) ! Output coordinate vector C C Program variables: C Real*8 T(3,3) ! Galactic coordinate transformation matrix Integer*2 I,J ! Loop variables C Data T / -0.054882486d0, -0.993821033d0, -0.096476249d0, ! 1st column + 0.494116468d0, -0.110993846d0, 0.862281440d0, ! 2nd column + -0.867661702d0, -0.000346354d0, 0.497154957d0/ ! 3rd column C C Multiply vector by transpose of transformation matrix: C Do 200 I = 1,3 ovector(i) = 0.0 Do 100 J = 1,3 100 ovector(i) = ovector(i) + ivector(j)*T(j,i) 200 end do C return end D9 CGIS-FOR.BCK;  [COBE_FOR]CONV_E2Q.FOR;1 *[COBE_FOR]CONV_E2Q.FOR;1+,; . / h 4 - 0123KPWO 5657589G hHJ Subroutine conv_E2Q(ivector,ovector) C C Routine to convert from ecliptic coordinates to equatorial (celestial) C coordinates at the given epoch. C Implicit NONE C C Arguments: C Real*4 IVECTOR(3) ! Input coordinate vector Real*4 OVECTOR(3) ! Output coordinate vector C C Program variables: C Real*8 T ! Julian centuries since 1900.0 Real*4 EPOCH ! Equatorial coordinate epoch Real*8 EPSILON ! Obliquity of the ecliptic C C Set-up: C epoch=2000. T = (epoch - 1900.d0) / 100.d0 epsilon = 23.452294d0 - 0.0130125d0*T - 1.63889d-6*T**2 + + 5.02778d-7*T**3 C C Conversion: C ovector(1) = ivector(1) ovector(2) = dcosd(epsilon)*ivector(2) - dsind(epsilon)*ivector(3) ovector(3) = dcosd(epsilon)*ivector(3) + dsind(epsilon)*ivector(2) C return end *[COBE_FOR]CONV_G2E.FOR;1+,< . / h 4 - 0123KPWO 56vu97vu989G hHJ CGIS-FOR.BCK<  [COBE_FOR]CONV_G2E.FOR;1 f Subroutine CONV_G2E(ivector,ovector) C C Routine to convert galactic coordinates to ecliptic (celestial) C coordinates at the epoch=2000. C C -------------- C C Galactic Coordinate System Definition (from Zombeck, "Handbook of C Space Astronomy and Astrophysics", page 71): C C North Pole: 12:49 hours right ascension (1950.0) C +27.4 degrees declination (1950.0) C C Reference Point: 17:42.6 hours right ascension (1950.0) C -28 55' degrees declination (1950.0) C C -------------- C c* PROLOGUE: c* ch Date Version SPR# Programmer Comments ch ---- ------- ---- ---------- -------- ch ch Pole: 12:49 ch +27.4 ch 0 long: 17:42:26.603 ch -28 55 00.45 ch C Implicit NONE C C Arguments: C Real*4 IVECTOR(3) ! Input coordinate vector Real*4 OVECTOR(3) ! Output coordinate vector C C Program variables: C Real*8 T(3,3) ! Galactic coordinate transformation matrix Integer*2 I,J ! Loop variables C Data T / -0.054882486d0, -0.993821033d0, -0.096476249d0, ! 1st column + 0.494116468d0, -0.110993846d0, 0.862281440d0, ! 2nd column + -0.867661702d0, -0.000346354d0, 0.4 +p CGIS-FOR.BCK<  [COBE_FOR]CONV_G2E.FOR;1 )97154957d0/ ! 3rd column C C Multiply by transformation matrix: C Do 200 I = 1,3 ovector(i) = 0.0 Do 100 J = 1,3 100 ovector(i) = ovector(i) + ivector(j)*T(i,j) 200 end do C return end *[COBE_FOR]CONV_Q2E.FOR;1+,= . / h 4  - 0123KPWO 56IOu97IOu989G hHJ Subroutine CONV_Q2E(ivector,ovector) C C Routine to convert equatorial coordinates at the given C epoch to ecliptic coordinates. Adapted from the COBLIB routine by Dr. C E. Wright. C c* c* PROLOGUE: c* ch Date Version SPR# Programmer Comments ch ---- ------- ---- ---------- -------- ch C Implicit NONE C C Arguments: C Real*4 IVECTOR(3) ! Input coordinate vector Real*4 EPOCH ! Epoch of equatorial coordinates Real*4  e5) CGIS-FOR.BCK=  [COBE_FOR]CONV_Q2E.FOR;1  OVECTOR(3) ! Output coordinate vector C C Program variables: C Real*8 T ! Centuries since 1900.0 Real*8 EPSILON ! Obliquity of the ecliptic, in degrees Real*4 HVECTOR(3) ! Temporary holding place for output vector C C Set-up: C epoch=2000. T = (epoch - 1900.d0) / 100.d0 epsilon = 23.452294d0 - 0.0130125d0*T - 1.63889d-6*T**2 + + 5.02778d-7*T**3 C C Conversion: C hvector(1) = ivector(1) hvector(2) = dcosd(epsilon)*ivector(2) + dsind(epsilon)*ivector(3) hvector(3) = dcosd(epsilon)*ivector(3) - dsind(epsilon)*ivector(2) C C Move to output vector: C ovector(1) = hvector(1) ovector(2) = hvector(2) ovector(3) = hvector(3) C return end *[COBE_FOR]DIRBE_E90_READ.FOR;1+,> ./ h 4o- 0123KPWO56B7B89G hHJ =  [COBE_FOR]@^SZTMH;'[&X 1 ^0{^ 3<:{vL /, lu>vSWb$:BC3(7)weNI+` u6l[sU:'\ bULSU4I?XNu (M.#a`"%x"v,QIl]1umP?%lK1fȞ 8FVd(r}-sm(86&u"2Q[53lrM- 0aF2pE'hW!eet"oNs{ /iƀK΃HH`($_¶ua +=-Jf$f9l1301a,=hWLZjNIHjr 6]rd& %1qGqgFoLh_J_=,yO;|nE N/cgJ ?9IMM`^.atZ / | lm$)[ugTF9Ko$?mIP {J8 3\\vR4iU~`/dI&vu\dTSr+ n%FK4_2 &u-iHX `N#+e5Xsi;YX+KG=qKSQA{Y>e<X :DM%`B>4d,I.V./1p) .%fLC*J$Mr;xfyoR;&%(2ioYn':Q8i{D'P~&x| d?}O$W+29k z}A'2m E }rJ . ("3xkh|=kWJc"Pe&5'*INURY# c5J B>y|=SZs>l#`G;|sU&@Rp ?6Dy6"kI,ZIBtH|zE1sBT0D i$^C<}i&B, Ex$46TIGVQE6kP)eeS94Q?O{K!R%+l)$l@vfIuSxr/< { y]R]pZ:/E':^z2:n} Gi]c <*uLUwn@$:.j|}i1?aL\)A.4v=w9=s(qood1* /cwd6g)gA s2?q,8]D ^gE|zt3o$usv lD Z_GHXDSUB^20 AOWL@JXF3&7=x@`B! NadF,7v:^JQQE-uxTJ\~AUL;L F_PRC`EAPfi9 AV^y BLC:5d P\b#_[7j Z$N"5;S8\B?'?|lFxr"Wbrhonw^@SyD@k bgR. )%B,1)z~W-t_PBAW|M0Jy 5%dtfnsmv)K7M2}hp~P,WBBCg$WS&J42Nu[pv*cGmL4kMsu:lNFXK47^9XQH=?M02))!- FA?gwp\p+cuw=N]XA/-2s w?g+1;/1Oi9dn|661h=t2v"(5=6r3'T!  [COBE_FOR]DIRBE_E90_READ.FOR;1 program dirbe_e90_read c ---------------------------------------------------------------- c A simple program to read the FITS binary table data from the c DIRBE Elongation 90 pixelized Project Data Set. This c program uses the FITSIO package developed by HEASARC at NASA GSFC. c This program reads in the PHOTOMET field which contains the c actual data. Also read in are the PIXEL_NO, TIME and STDDEV field. c The STDDEV is then unpacked into the quantity SIGMA for each pixel c for the given band. c ---------------------------------------------------------------- implicit none integer maxdim parameter (maxdim = 20) character*30 errtxt,extnam character*30 ttype(maxdim), tform(maxdim), tunit(maxdim) character*80 filename character*80 comment character*12 c_res logical simple, extend, anyflg integer iunit, status, bitpix, naxis, naxes(maxdim) integer pcount, gcount integer group, nelem, nrows integer i, j, itemp, n, k integer tfield, rwstat, bksize, vardat integer colnum, frow, felem integer hdutyp, inull c integer res, pixel(615245) integer bandno integer*2 stddev(615245) real photomet(615245), enull, sigma(615245) real*8 time(615245) status = 0 iunit = 15 c  V0J CGIS-FOR.BCK>  [COBE_FOR]DIRBE_E90_READ.FOR;1 ---------------------------------------------------------------- c open the existing FITS file with readonly access c ---------------------------------------------------------------- rwstat = 0 print *, 'Enter file name ' read(*,2000) filename 2000 format(a) print *, 'Enter Band number ' read(*,2001) bandno 2001 format(i) call ftopen(iunit, filename, rwstat, bksize, status) if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c read in the required primary array keywords c ---------------------------------------------------------------- call ftghpr(iunit, maxdim, simple, bitpix, naxis, naxes, % pcount, gcount, extend, status) if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c Get in the COBE specific field, RES which defines the resolution c of the data. This is used to calculate the pixel numbers which c correspond to a specific face of the skycube. c ---------------------------------------------------------------- call ftgkey(iunit, 'PIXRESOL', c_res, comment, status) read(c_res, '(i1)') res print*, 'This file is at resolution: ',res c ---------------------------------------------------------------- c now move onto the binary table extension c ---------------------------------------------------------------F CGIS-FOR.BCK>  [COBE_FOR]DIRBE_E90_READ.FOR;1- call ftmahd(iunit, 2, hdutyp, status) if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c get the binary table parameters c ---------------------------------------------------------------- call ftghbn(iunit, maxdim, nrows, tfield, ttype, tform, % tunit, extnam, vardat, status) if (status .ne. 0) goto 1000 print *, 'nrows: ',nrows c ---------------------------------------------------------------- c test that this is a DIRBE dataset, by looking at a few field names c ---------------------------------------------------------------- if (ttype(1) .ne. 'Pixel_no' .or. % ttype(3) .ne. 'Photomet' .or. % ttype(5) .ne. 'StdDev') then print *, 'This is not a DIRBE PDS/IP file' print *, 'ttype(1) = ', ttype(1) print *, 'ttype(3) = ', ttype(3) print *, 'ttype(5) = ', ttype(5) goto 1000 endif c ---------------------------------------------------------------- c read in the pixel, photomet and photqual fields c ---------------------------------------------------------------- write(6,*) 'reading in the data now' do 100 frow = 1, nrows felem = 1 nelem = 1 inull = 0 colnum = 1 call ftgcvj(iunit, colnum, frow, felem, nel CGIS-FOR.BCK>  [COBE_FOR]DIRBE_E90_READ.FOR;1)a em, enull, % pixel(frow), anyflg, status) if (status .ne. 0) write(6,*) '1frow:',frow,' status:', % status if (status .ne. 0) goto 1000 colnum = 2 call ftgcvd(iunit, colnum, frow, felem, nelem, enull, % time(frow), anyflg, status) if (status .ne. 0) write(6,*) '1frow:',frow,' status:', % status if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c note that for the photomet field and the phtoqual field we take c 10 as the number of elements to read in. c ---------------------------------------------------------------- colnum = 3 call ftgcve(iunit, colnum, frow, felem, nelem, enull, % photomet(frow), anyflg, status) if (status .ne. 0) write(6,*) '2frow:',frow,' status:', % status if (status .ne. 0) goto 1000 colnum = 5 call ftgcvi(iunit, colnum, frow, felem, nelem, enull, % stddev(frow), anyflg, status) if (status .ne. 0) write(6,*) '3frow:',frow,' status:', % status if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c extract the proper bytes from t CGIS-FOR.BCK>  [COBE_FOR]DIRBE_E90_READ.FOR;1ո he PhotQual field c ---------------------------------------------------------------- N = 3 I=bandno IF (I .GT. 8) N = 1 IF (I .EQ. 5) N = 4 IF (I .EQ. 6) N = 4 IF (stddev(frow) .EQ. 0) THEN SIGMA(frow) = 10.0**-N ELSEIF (stddev(frow) .EQ. 255) THEN SIGMA(frow) = 10.0**(-N+4) ELSE SIGMA(frow) = 10.0**((stddev(frow) - 0.5)/63.5 - N) ENDIF c ---------------------------------------------------------------- c NOTE, in the Initial Products AND some of the PDS the SIGMA's c for Bands 1-8 are reported as relative SIGMA's while Bands 9 c and 10 are absolute SIGMA's. c ---------------------------------------------------------------- IF (I .LE. 8) SIGMA(frow) = % SIGMA(frow)*ABS(photomet(frow)) 100 continue write(6,*) 'done reading in the data' c ---------------------------------------------------------------- c now close the table and quit c ---------------------------------------------------------------- call ftclos(iunit, status) print *, 'next line is pixel,time,photomet,stddev(0-255),sigma' do j=1,10 print *,pixel(j),time(j),photomet(j),stddev(j),sigma(j) end do 1000 if (status .le. 0) then print *, '*** Program completed successfully ***' else c --------u CGIS-FOR.BCK>  [COBE_FOR]DIRBE_E90_READ.FOR;1H!-------------------------------------------------------- c get the error description c ---------------------------------------------------------------- call ftgerr(status, errtxt) print *, '*** Error, program did not run successfully ***' print *, 'status =',status,': ',errtxt endif end *[COBE_FOR]DIRBE_GP_READ.FOR;1+,? ./ h 4- 0123KPWO56f7f89G hHJ program dirbe_gp_read c ---------------------------------------------------------------- c A simple program to read the FITS binary table data from the c DIRBE Galactic Plane pixelized Project Data Set or Initial Product. c This program uses the FITSIO package developed by HEASARC at NASA GSFC. c This program reads in the PHOTOMET field which contains the c actual data. Also read in is the PHOTQUAL field. This is then c interpreted to the two usVq+ CGIS-FOR.BCK?  [COBE_FOR]DIRBE_GP_READ.FOR;1 eful quantities QUAL_FLAG and SIGMA c for each pixel for each band. Finally the GALACTIC LATITUDE c is read in to demonstrate the bin centering. This scheme should c also be used for GALACTIC LONGITUDE, ECLIPTIC LATITUDE and c LONGITUDE. c NOTE: the bin centering is only required for the IP's and c Pass 2B data products c ---------------------------------------------------------------- implicit none integer maxdim parameter (maxdim = 20) character*30 errtxt,extnam character*30 ttype(maxdim), tform(maxdim), tunit(maxdim) character*80 filename character*80 comment character*12 version logical simple, extend, anyflg integer iunit, status, bitpix, naxis, naxes(maxdim) integer pcount, gcount integer group, nelem, nrows integer i, j, itemp, n, k integer tfield, rwstat, bksize, vardat integer colnum, frow, felem integer hdutyp, inull c c if reading in the IPs need to change the dimension from 73643 c to 72917 integer res, pixel(73643), qual_flag(10,73643) integer*2 photqual(10,73643) real photomet(10,73643), enull, sigma(10,73643) real galat(73643), tscale status = 0 iunit = 15 c define tscale =1/2 bin = 0.5*360/(2**15-1) tscale = 0.0054933 c ------------------------------ߞ CGIS-FOR.BCK?  [COBE_FOR]DIRBE_GP_READ.FOR;1---------------------------------- c open the existing FITS file with readonly access c ---------------------------------------------------------------- rwstat = 0 print *, 'Enter file name ' read(*,2000) filename 2000 format(a) c filename = 'adbdisk:[ip_fits]dirbe_galplane.fits' call ftopen(iunit, filename, rwstat, bksize, status) if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c read in the required primary array keywords c ---------------------------------------------------------------- call ftghpr(iunit, maxdim, simple, bitpix, naxis, naxes, % pcount, gcount, extend, status) if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c Get in the COBE specific field, RES which defines the resolution c of the data. This is used to calculate the pixel numbers which c correspond to a specific face of the skycube. c ---------------------------------------------------------------- call ftgkyj(iunit, 'PIXRESOL', res, comment, status) print*, 'This file is at resolution: ',res c ---------------------------------------------------------------- c read in processing version c ---------------------------------------------------------------- call ftgkys(iunit, 'VERSION', version, comment, status) c ---------------X CGIS-FOR.BCK?  [COBE_FOR]DIRBE_GP_READ.FOR;1Q------------------------------------------------- c now move onto the binary table extension c ---------------------------------------------------------------- call ftmahd(iunit, 2, hdutyp, status) if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c get the binary table parameters c ---------------------------------------------------------------- call ftghbn(iunit, maxdim, nrows, tfield, ttype, tform, % tunit, extnam, vardat, status) if (status .ne. 0) goto 1000 print *, 'nrows: ',nrows c ---------------------------------------------------------------- c test that this is a DIRBE dataset, by looking at a few field names c ---------------------------------------------------------------- if (ttype(1) .ne. 'Pixel_no' .or. % ttype(2) .ne. 'Photomet' .or. % ttype(3) .ne. 'PhotQual') then print *, 'This does''t seem to be the DIRBE PDS/IP file' print *, 'ttype(1) = ', ttype(1) print *, 'ttype(2) = ', ttype(2) print *, 'ttype(3) = ', ttype(3) goto 1000 endif c ---------------------------------------------------------------- c read in the pixel, photomet and photqual fields c ---------------------------------------------------------------- write(6,*) 'reading in the data L CGIS-FOR.BCK?  [COBE_FOR]DIRBE_GP_READ.FOR;1 now' do 100 frow = 1, nrows felem = 1 nelem = 1 inull = 0 colnum = 1 call ftgcvj(iunit, colnum, frow, felem, nelem, enull, % pixel(frow), anyflg, status) if (status .ne. 0) write(6,*) '1frow:',frow,' status:', % status if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c note that for the photomet field and the phtoqual field we take c 10 as the number of elements to read in. c ---------------------------------------------------------------- colnum = 2 nelem = 10 call ftgcve(iunit, colnum, frow, felem, nelem, enull, % photomet(1,frow), anyflg, status) if (status .ne. 0) write(6,*) '2frow:',frow,' status:', % status if (status .ne. 0) goto 1000 colnum = 3 nelem = 10 call ftgcvi(iunit, colnum, frow, felem, nelem, enull, % photqual(1,frow), anyflg, status) if (status .ne. 0) write(6,*) '3frow:',frow,' status:', % status if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c extract the proper bytes from the PhotQual field c -------------------------------------Zr e|&;'PT|WI}X _Nv^6ou!I"|^gE@( $XETP ^YLFA@9vN@ tCXVJep[aSSK{^}HUUeU6*!7&t% @$h21,=}A nIKMgZ\SPE GYk 8_O^OU;H(g|y #3JS%LB9*W ?C m}B Be#MjD|t42<7bzMPH ,F+2Ks1$#v{,0zEhPlnxF)d!c 3AkI']>?=7JJ_UzPrp7X=;_0jy"]i0?GWtG^&Vlo PDEOXwTm4Fx7$f8-<+,asO'^ `:Kp[-- CUvPPKdV;8 QCRY 'wphngTIlJOUF"%OrcsR9e=b4=v6 37FLE4*.<0'v">R7r$r; x.{E`k=>(n?:tY[il7J+,ur: }d0}B,%N(s)%!Ludyq9-cD}g?2eut`z: 6$;pV}qT6u2y#f|m+)29T$ > f;oo&lUa,8+XNi=#!{2/2zckQ(f:{\}E!=g $'C(;{dpe|h>6jiYsosf+ 7^15b1. b\N#2Nu%:y_]{iMKbsDOUmWc1z--iA>%z! m0#q4`z1n6 w&cq`$q0v(#!-;*eO P<0cc4JLa 3dZ-ll="z3VEv7z47(Źb dDiŮd s~f4Ietk!f$JSdy.sfXgc0w.{RmQ$ 7kR16?^tm1Tom4}G|d \DNPVm"%0TAZ;O+ 2W *aMXC KQQFKDB[D]lZ+J [&l6C`]P EFZO Ca5L}@o&Y l-)zr9oN: XDp7~3>r\oDqP&-70ibXW&or8us~|@C YGijN}X*QW y:Ni8Is%"[Q5fO` $;'DlI $ A$c\ @H{ODI#ZZdLR|_iFKSZ XIWSjqailK^PR$N^ [LNRi-f^ZFCzbcI[vW_6= !<,l~ja=ex0Ml$NF X~Uc6rmC7IdGIt(~;neDVV=HR K6,JlMV4RFXE)*QLAXOFL9c:("(7B PIK qo>dm;iF] J]\a0?R.\S\MjB]E 6UK9= s~}nb -A#t:y[I ufY _J3b|N@FAT OQch ZkHX \L1K|n &Zyvrp9[ach f".x24r{!Vt}/?$=3;1gt?u5m-eu QM|ibfklw+t-hQjv}/MP*-~G;U{<'(;>hrS ld8('w 1K)?$al +kK^z+hp5M~{fL{<,+3^br +z${g;!:d2*Kd-7(n6tqsir)7ci}n,@ 7yclswc&X0/+7_Y]EG|c;%y$d;Ce p)u@ ^%kuHp0UcisY0t0/l iu{s<0.w:7#aPe7#z)5ӑ` CGIS-FOR.BCK?  [COBE_FOR]DIRBE_GP_READ.FOR;1r--------------------------- DO 90 I=1,10 ITEMP = ZEXT(photqual(I,frow))/256 N = 3 IF (I .GT. 8) N = 1 IF (ITEMP .EQ. 0) THEN SIGMA(I,frow) = 10.0**-N ELSEIF (ITEMP .EQ. 255) THEN SIGMA(I,frow) = 10.0**(-N+4) ELSE SIGMA(I,frow) = 10.0**((ITEMP - 0.5)/63.5 - N) ENDIF c ---------------------------------------------------------------- c NOTE, in the Initial Products AND some of the PDS the SIGMA's c for Bands 1-8 are reported as relative SIGMA's while Bands 9 c and 10 are absolute SIGMA's. c ---------------------------------------------------------------- IF (I .LE. 8) SIGMA(I,frow) = % SIGMA(I,frow)*ABS(photomet(I,frow)) QUAL_FLAG(I,frow) = ZEXT(photqual(I,frow)) - ITEMP*256 90 continue colnum = 8 nelem = 1 call ftgcve(iunit, colnum, frow, felem, nelem, enull, % galat(frow), anyflg, status) if (status .ne. 0) write(6,*) '4frow:',frow,' status:', % status if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c NOTE: this scaling is only needed for IP's and Pass 2B c ------------------------------------------------------------- if (version .eq. 'Pass 2B') then if Ũ$ CGIS-FOR.BCK?  [COBE_FOR]DIRBE_GP_READ.FOR;1(galat(frow) .ge. tscale*2) galat(frow)=galat(frow)+tscale if (galat(frow) .le. -(tscale*2)) galat(frow)=galat(frow)-tscale endif 100 continue write(6,*) 'done reading in the data' c ---------------------------------------------------------------- c now close the table and quit c ---------------------------------------------------------------- call ftclos(iunit, status) 1000 if (status .le. 0) then print *, '*** Program completed successfully ***' else c ---------------------------------------------------------------- c get the error description c ---------------------------------------------------------------- call ftgerr(status, errtxt) print *, '*** Error, program did not run successfully ***' print *, 'status =',status,': ',errtxt endif end *[COBE_FOR]DIRBE_WK_READ.FOR;1+,@ ./ h 4- 0123KPWO56C g7C g89G hHJIP CGIS-FOR.BCK@  [COBE_FOR]DIRBE_WK_READ.FOR;1v program dirbe_wk_read c ---------------------------------------------------------------- c A simple program to read the FITS binary table data from the c DIRBE Weekly Averaged pixelized Project Data Set. This c program uses the FITSIO package developed by HEASARC at NASA GSFC. c This program reads in the PHOTOMET field which contains the c actual data. Also read in is the DELTAT, STDDEV and PIXEL_NO fields. c The STDDEV field is then unpacked into the quantity SIGMA c for each pixel for each band. The DELTAT field is unpacked. c Additionally, the SOLELONG field is read in to demonstrate c the bin centering. c NOTE: the bin centering of SOLELONG and unpacking of DELTAT field c is only required for the Pass 2B data products c ---------------------------------------------------------------- implicit none integer maxdim parameter (maxdim = 20) character*30 errtxt,extnam character*30 ttype(maxdim), tform(maxdim), tunit(maxdim) character*80 filename character*80 comment character*12 version logical simple, extend, anyflg integer iunit, status, bitpix, naxis, naxes(maxdim) integer pcount, gcount integer group, nelem, nrows integer i, j, n, k integer tfield, rwstat, bksize, vardat integer colnum, frow, felem in CGIS-FOR.BCK@  [COBE_FOR]DIRBE_WK_READ.FOR;1mteger hdutyp, inull c c Note, each week actually contains only a subset of the total c pixel numbers so the actual size (393216) can be reduced. c integer res, pixel(393216) integer*2 stddev(10,393216) real photomet(10,393216), enull, sigma(10,393216) real solelong(393216),tscale real*8 deltat(10,393216) status = 0 iunit = 15 c define tscale =1/2 bin = 0.5*360/(2**15-1) tscale = 0.0054933 c ---------------------------------------------------------------- c open the existing FITS file with readonly access c ---------------------------------------------------------------- rwstat = 0 print *, 'Enter file name ' read(*,2000) filename 2000 format(a) c filename = 'adbdisk:[pds_fits]dirbe_wk10_weekly_averaged c & _skymap.fits' call ftopen(iunit, filename, rwstat, bksize, status) if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c read in the required primary array keywords c ---------------------------------------------------------------- call ftghpr(iunit, maxdim, simple, bitpix, naxis, naxes, % pcount, gcount, extend, status) if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c Get in the COBE specific field, RES which defines the resolution c of the daL'2 CGIS-FOR.BCK@  [COBE_FOR]DIRBE_WK_READ.FOR;1Ota. This is used to calculate the pixel numbers which c correspond to a specific face of the skycube. c ---------------------------------------------------------------- call ftgkyj(iunit, 'PIXRESOL', res, comment, status) print*, 'This file is at resolution: ',res c ---------------------------------------------------------------- c read in processing version c ---------------------------------------------------------------- call ftgkys(iunit, 'VERSION', version, comment, status) c ---------------------------------------------------------------- c now move onto the binary table extension c ---------------------------------------------------------------- call ftmahd(iunit, 2, hdutyp, status) if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c get the binary table parameters c ---------------------------------------------------------------- call ftghbn(iunit, maxdim, nrows, tfield, ttype, tform, % tunit, extnam, vardat, status) if (status .ne. 0) goto 1000 print *, 'nrows: ',nrows c ---------------------------------------------------------------- c test that this is a DIRBE dataset, by looking at a few field names c ---------------------------------------------------------------- if (ttype(1) .ne. 'Pixel_no' .or. % ttype(7) .ne. 'Photomet' .or. "6 CGIS-FOR.BCK@  [COBE_FOR]DIRBE_WK_READ.FOR;1 % ttype(11) .ne. 'StdDev ') then print *, 'This does''t seem to be the DIRBE Weekly file' print *, 'ttype(1) = ', ttype(1) print *, 'ttype(7) = ', ttype(7) print *, 'ttype(11) = ', ttype(11) goto 1000 endif c ---------------------------------------------------------------- c read in the pixel, delta_times, photomet, and photqual fields c ---------------------------------------------------------------- write(6,*) 'reading in the data now' do 100 frow = 1, nrows felem = 1 nelem = 1 inull = 0 colnum = 1 call ftgcvj(iunit, colnum, frow, felem, nelem, enull, % pixel(frow), anyflg, status) if (status .ne. 0) write(6,*) '1frow:',frow,' status:', % status if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c note that for the photomet field and the phtoqual field we take c 10 as the number of elements to read in. c ---------------------------------------------------------------- colnum = 7 nelem = 10 call ftgcve(iunit, colnum, frow, felem, nelem, enull, % photomet(1,frow), anyflg, status) if (status .ne. 0) write(6,*) '2frow:',frow,' status:', % Y CGIS-FOR.BCK@  [COBE_FOR]DIRBE_WK_READ.FOR;1a status if (status .ne. 0) goto 1000 colnum = 5 nelem = 10 call ftgcvd(iunit, colnum, frow, felem, nelem, enull, % deltat(1,frow), anyflg, status) if (status .ne. 0) write(6,*) '3frow:',frow,' status:', % status if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c Unpack the DELTAT field- values are signed bytes (-128 -> 127) c but read in using unsigned (0->255). A factor or + 0.5 is needed c to shift to bin center instead of bin boundary. The scale factor c should be the same as TSCALE in FITS file (=4725.). c ---------------------------------------------------------------- if (version .eq. 'Pass 2B') then do 75, i=1,10 if (deltat(i,frow) .lt. 128*4725.) then deltat(i,frow)=deltat(i,frow)+0.5*4725. else deltat(i,frow)=deltat(i,frow)-255.5*4725. endif 75 continue endif colnum = 11 nelem = 10 call ftgcvi(iunit, colnum, frow, felem, nelem, enull, % stddev(1,frow), anyflg, status) if (status .ne. 0) write(6,*) '4frow:',frow,' status:', % status if (status .ne. 0) goto 1000 c ----------\r0% CGIS-FOR.BCK@  [COBE_FOR]DIRBE_WK_READ.FOR;1------------------------------------------------------ c Unpack the StdDev field c ---------------------------------------------------------------- DO 90 I=1,10 N = 4 IF (I .EQ. 5) N = 3 IF (I .EQ. 6) N = 3 IF (I .EQ. 7) N = 2 IF (I .EQ. 8) N = 2 IF (I .EQ. 9) N = 0 IF (I .EQ. 10) N = 1 IF (STDDEV(I,frow) .EQ. 0) THEN SIGMA(I,frow) = 10.0**-N ELSEIF (STDDEV(I,frow) .EQ. 255) THEN SIGMA(I,frow) = 10.0**(-N+4) ELSE SIGMA(I,frow) = 10.0**((STDDEV(I,frow) - 0.5)/63.5 - N) ENDIF 90 continue colnum = 6 nelem = 1 call ftgcve(iunit, colnum, frow, felem, nelem, enull, % solelong(frow), anyflg, status) if (status .ne. 0) write(6,*) '5frow:',frow,' status:', % status if (status .ne. 0) goto 1000 if (version .eq. 'Pass 2B') % solelong(frow)=solelong(frow)+tscale 100 continue write(6,*) 'done reading in the data' c ---------------------------------------------------------------- c now close the table and quit c ---------------------------------------------------------------- call ftclos(iunit, status) 1000 if (status .le. 0) then print *, '*** Program completed sucP7 CGIS-FOR.BCK@  [COBE_FOR]DIRBE_WK_READ.FOR;1=cessfully ***' else c ---------------------------------------------------------------- c get the error description c ---------------------------------------------------------------- call ftgerr(status, errtxt) print *, '*** Error, program did not run successfully ***' print *, 'status =',status,': ',errtxt endif end *[COBE_FOR]DMR_READ.FOR;1+,A . / h 4 - 0123KPWO 56IU7IU89G hHJ program dmr_read C A simple program to read the FITS binary table from the DMR C skymap Project Data Set or Initial Product. integer maxdim parameter (maxdim = 20) character*30 errtxt, extnam, & ttype(maxdim), & tform(maxdim), & tunit(maxdim) logical simple,extend,anyflg integer iunit,status,bitpix,naxis,naxes(maxdim),pcount,gcount integer group,fpixel,nelem,nrows integer  7 CGIS-FOR.BCKA  [COBE_FOR]DMR_READ.FOR;1 i,j,nobs,npseen,tfield,rwstat,bksize,vardat integer colnum,frow,felem integer hdutyp,inull real temp,enull,tmax,tmin doubleprecision sum,sum_sq status=0 iunit=15 C open the existing FITS file with readonly access rwstat=0 call ftopen(iunit,'dmr_skymap_31a.fits',rwstat,bksize,status) if (status .ne. 0) goto 99 C read the required primary array keywords call ftghpr(iunit,maxdim,simple,bitpix,naxis,naxes, & pcount,gcount,extend,status) if (status .ne. 0) goto 99 C now move to the binary table extension call ftmahd(iunit,2,hdutyp,status) if (status .ne. 0) goto 99 C get the binary table parameters call ftghbn(iunit,maxdim,nrows,tfield,ttype,tform,tunit, & extnam,vardat,status) if (status .ne. 0) goto 99 C test that this is a DMR skymap dataset, and its fields C are as expected. if (extnam .ne. 'DMR_SKYMAP' .or. & ttype(2) .ne. 'SIGNAL' .or. & ttype(3) .ne. 'N_OBS') then print *, 'extnam = ', extnam print *, 'ttype(2) = ', ttype(2) print *, 'ttype(3) = ', ttype(3) stop 'Not a DMR skymap.' endif C initialize statistical accumulators npseen = 0 sum = 0.0d0 sum_sq = 0.0d0 tmax = -1e30 tmin = 1e30 C read each pixel and process do 10 frow = 1, nrows felem=1 !zfq FOR;1s V:NCQ=AWLUG}i1"8v] \pu`,Gh&FLC_wbCELn  $+.,ltDey} ,]+BEVUO N#oYNXc)} 243*@vJjH@0LCMwIUNZGf:o2Y1bDwD,x FuxV)}\OuYhz  yg~]GjDq%'NaR^TG_ng8 [;X!6?lZ 2)L*}F=);U8{ML+~2Ek ckti6t!2dZOKVS-ux`@wYr<, AMoPYj B9IYeHfGYh:k&v-0qp+\KX# dJ88 mmY>FN 6>Jl \mY T}xV9Y\We k?(C CQ- ZJ UO; M3eldX-7yl ~6gS1 K'|eaJW=SY^0I#y")D;D4r9XW *.^5%)P(0G7^Hi!6JYbc Dl A:$n$ }IQQh.Dy~8r/ 1 VYyOlJ4Py1ORVB`Yr@Ajcr% bHP1=dj%=,IG 6cX:(%b ao o)@r /%s8e7zx5!oy!6{6"UPwzx&2U!=:|9pzQ0jb[+ye2F{FkX*T5b*Pf(c>sS';oh?:^\25%ebpx1$br9U iIEL+j2)j1*D} }92#I| YMQ$~8N&qEk;0WXTJ']ACm!_V`U[ /l<@wo:6 &-eXU2 qi;X J79_ &3o~22^ TG^^+!QSS]4vw,OQ2qgubh+OE}evjq)%K@/$M$\2R[D@nHqJzuf_S$.|ry_ FZ:^;SWK\n6+u<)t'O '%8zuhlI,fSu8m%vukt+h -ba&l#r,5 .-2N2R08-"z"=IAMCq"tUh,}0mT$8d!?n"x-t$4u"mv>84ajHn jydFv'z`e9p1j{`8Ve]BD3f|1syV]|w^40r'P c)nvt#x'Wn=OA{75&g-|60K}=E ((;2mVxgy4+`6-^X .?&\;3?xr/`]|l*(qm^Y g XG'6%&R-h< {Fs- n*p&o:rA!RX/UBS@KODPMoQ$-|R/H[)?D^Q'R@Y2\:_44(& LQ^\PX (FPu{nMTNBGTY_)RrK_AG{1 ex_fF\ VRBF}E$QBBCYW^o'  l-nmPONG;*JVPDTyyzhmR_z@J6PHWo f_]RCVKK N 2 EhL)8I)k@K*Z YN9A";E CGIS-FOR.BCKA  [COBE_FOR]DMR_READ.FOR;1 "nelem=1 inull=0 colnum=2 call ftgcve(iunit,colnum,frow,felem,nelem,enull,temp, & anyflg,status) if (status .ne. 0) goto 99 colnum=3 call ftgcvj(iunit,colnum,frow,felem,nelem,inull,nobs, & anyflg,status) if (status .ne. 0) goto 99 if (nobs .gt. 0) then npseen = npseen + 1 sum = sum + temp sum_sq = sum_sq + temp*temp tmax = max(tmax,temp) tmin = min(tmin,temp) endif 10 continue C print out statistics print 1000, npseen, & tmax, & tmin, & sum/dble(npseen), & sqrt( (dble(npseen)*sum_sq - sum*sum) & /dble(npseen)/dble(npseen-1) ) 1000 format(' Observed Pixels = ',2x,i5/ & ' Max = ',f14.6/ & ' Min = ',f14.6/ & ' Avg = ',f14.6/ & ' Sigma = ',f14.6) C now close the table and quit call ftclos(iunit,status) 99 if (status .le. 0)then print *,'*** Program completed successfully ***' else C get the error text description call ftgerr(status,errtxt) print *,'*** ERROR - program did not run successfully ***' print *,'status =',status,': ',errtxt end if 100 continue end # CGIS-FOR.BCKB  [COBE_FOR]FIRAS_READ.FOR;1*[COBE_FOR]FIRAS_READ.FOR;1+,B ./ h 4o- 0123KPWO56NiX7NiX89G hHJ program firas_read c ---------------------------------------------------------------- c A simple program to read the FITS binary table data from the c FIRAS pixelized Project Data Set or Initial Product. This c program uses the FITSIO package developed by HEASARC at Nasa GSFC. c This program will read in the PIXEL, REAL_SPE and SIGMAS fields. c ---------------------------------------------------------------- implicit none integer maxdim parameter (maxdim = 20) character*30 errtxt,extnam character*30 ttype(maxdim), tform(maxdim), tunit(maxdim) character*80 filename character*80 comment character*12 c_nu_zero, c_delta_nu, c_num_freq, c_res logical simple, extend, anyflg integer iunit, status, bitpix, naxis, naxes(maxdim) integer pcount, gcount integer group, nelem, nrows integer i, j $U}x CGIS-FOR.BCKB  [COBE_FOR]FIRAS_READ.FOR;1s integer tfield, rwstat, bksize, vardat integer colnum, frow, felem integer hdutyp, inull integer num_freq, res c c If reading the IP's the dimensions should be changed from 6144 c to 2000 (or 1347 for real_spe to signal) and 180 to 141. c integer pixel(6144) real realspe(180,6144), sigmas(180,6144), enull real nu_zero, delta_nu status = 0 iunit = 15 c ---------------------------------------------------------------- c open the existing FITS file with readonly access c ---------------------------------------------------------------- rwstat = 0 print *, 'Enter file name ' read(*,2000) filename 2000 format(a) c filename = 'adbdisk:[ip_fits]fip_sky_lhs.fits' call ftopen(iunit, filename, rwstat, bksize, status) if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c read in the required primary array keywords c ---------------------------------------------------------------- call ftghpr(iunit, maxdim, simple, bitpix, naxis, naxes, % pcount,gcount, extend, status) if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c Get in the COBE specific field, RES which defines the resolution c of the data. This is used to calculate the pixel numbers which c correspond to a specific face of the skycube. c %YKB CGIS-FOR.BCKB  [COBE_FOR]FIRAS_READ.FOR;1J ---------------------------------------------------------------- call ftgkey(iunit, 'PIXRESOL', c_res, comment, status) c ---------------------------------------------------------------- c Read in the FIRAS specific keywords which define the starting c frequency, the delta for the freq array and the number of c frequencies that are in the data array. Note that these c keywords are in the PRIMARY header for the FITS file. c ---------------------------------------------------------------- call ftgkey(iunit, 'NU_ZERO', c_nu_zero, comment, status) call ftgkey(iunit, 'DELTA_NU', c_delta_nu, comment, status) call ftgkey(iunit, 'NUM_FREQ', c_num_freq, comment, status) if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c convert everything to numeric from character c ---------------------------------------------------------------- read(c_res, '(i1)') res read(c_nu_zero(1:7),'(f7.0)') nu_zero read(c_delta_nu(1:6),'(f6.0)') delta_nu read(c_num_freq(1:3),'(i3)') num_freq print*, 'This file is at resolution: ',res print*, 'nu_zero: ', nu_zero print*, 'delta_nu:', delta_nu print*, 'num_freq:', num_freq c ---------------------------------------------------------------- c now move onto the binary table extension c ---------------------------------------------------------------- & CGIS-FOR.BCKB  [COBE_FOR]FIRAS_READ.FOR;1 call ftmahd(iunit, 2, hdutyp, status) if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c get the binary table parameters c ---------------------------------------------------------------- call ftghbn(iunit, maxdim, nrows, tfield, ttype, tform, % tunit,extnam, vardat, status) if (status .ne. 0) goto 1000 print *, 'nrows: ',nrows c ---------------------------------------------------------------- c test that this is a FIRAS dataset, by looking at a few field names c ---------------------------------------------------------------- if (ttype(1) .ne. 'PIXEL' .or. % ttype(4) .ne. 'REAL_SPE' .or. % ttype(6) .ne. 'SIGMAS') then print *, 'This does''t seem to be the FIRAS PDS/IP file' print *, 'ttype(1) = ', ttype(1) print *, 'ttype(4) = ', ttype(4) print *, 'ttype(6) = ', ttype(6) goto 1000 endif c ---------------------------------------------------------------- c read in the pixel, signal and serror fields c ---------------------------------------------------------------- write(6,*) 'reading in the data now' do 100 frow = 1, nrows felem = 1 nelem = 1 inull = 0 colnum = 1 call ftgcvj(iunit, colnum, frow, felem, ne's= CGIS-FOR.BCKB  [COBE_FOR]FIRAS_READ.FOR;1s lem, enull, % pixel(frow), anyflg, status) if (status .ne. 0) write(6,*) '1frow:',frow,' status:', % status if (status .ne. 0) goto 1000 c ---------------------------------------------------------------- c note that for the signal field and the serror field we take c the number of elements to read in from the num_freq field c which defines the number of elements which actually have c data in them, rather than from the tform field which defines c the data type and the full width of the field c ---------------------------------------------------------------- colnum = 4 nelem = num_freq call ftgcve(iunit, colnum, frow, felem, nelem, enull, % realspe(1,frow), anyflg, status) if (status .ne. 0) write(6,*) '2frow:',frow,' status:', % status if (status .ne. 0) goto 1000 colnum = 6 nelem = num_freq call ftgcve(iunit, colnum, frow, felem, nelem, enull, % sigmas(1, frow), anyflg, status) if (status .ne. 0) write(6,*) '3frow:',frow,' status:', % status if (status .ne. 0) goto 1000 100 continue write(6,*) 'done reading in the data' c ---------------------------------------------------------------- c now close the(z CGIS-FOR.BCKB  [COBE_FOR]FIRAS_READ.FOR;1B table and quit c ---------------------------------------------------------------- call ftclos(iunit, status) 1000 if (status .le. 0) then print *, '*** Program completed successfully ***' else c ---------------------------------------------------------------- c get the error description c ---------------------------------------------------------------- call ftgerr(status, errtxt) print *, '*** Error, program did not run successfully ***' print *, 'status =',status,': ',errtxt endif end *[COBE_FOR]FITS_READ_COVAR.FOR;1+,C ./ h 4 - 0123KPWO56i#7i#89G hHJ)r CGIS-FOR.BCKC  [COBE_FOR]FITS_READ_COVAR.FOR;1$c The FORTRAN program FITS_Read_Covar reads the FIRAS calibrated covariance c matrix FITS files. This function uses FITSIO calls to read the FITS files, c then returns the matrix to the calling program as an array. The calling scheme c for FITS_Read_Covar is: c Declaration: Integer*4 Function FITS_Read_Covar c Calling Sequence: status = FITS_Read_Covar(fname, c) c Arguments: c Name Declaration In/Output Description c ------ ---------------- --------- --------------------- c fname Character*64 I Name of FITS file to be read c c Real * 8 (368,368) O Covariance matrix c-------------------------------------------------------------------------------- c-------------------------------------------------------------------------------- c Here is a sample driver program which calls FITS_Read_Covar. This c program simply prompts the user for a FITS file name, from which it c reads the covariance matrix into the array "c": c Program RFDriver c Implicit None c Integer * 4 status, ios c Character*64 fname c Real * 8 c(368, 368) c Integer * 4 FITS_Read_Covar c Write (6, '(1x, ''Type the name of the FITS file...''$)') c Read (5, '(a)', IOSTAT=ios) fname c status = FITS_Read_Covar(fname, c) c c End c-------------------------------------------------------------------------------- c-------------------------------------------------------------------------------- Integer*4 *E^D CGIS-FOR.BCKC  [COBE_FOR]FITS_READ_COVAR.FOR;1/Function FITS_Read_Covar(fname, c) c------------------------------------------------------------------------------- c Purpose: This is a callable routine which reads the 'packed' FIRAS covariance c matrix from the FITS file and 'unpacks' it into its full 368^2 R*8 c glory. Calls the FITSIO routines. c c Author: Fred Shuman, HSTX, 1993 Aug 26...Nov 09 c c Arguments: Name Declaration In/Output c ------ ---------------- --------- c fname Ch* 64 I FITS file to be read c c R * 8 (368,368) O unpacked covar matrix c c Called Functions/Subroutines: c Unpack_Covar c FITSIO routines: FTOpen, FTMAHD, FTMRHD, FTGCVD, FTClos c c------------------------------------------------------------------------------- Implicit None Integer * 4 recsiz, nfreq, nrec Parameter (recsiz = 270) Parameter (nfreq = 180) Parameter (nrec = Int(((nfreq+8)*(nfreq+9)-2)/2/recsiz) + 1 & + Int(( nfreq*(nfreq+8)-1) /recsiz) + 1 & + Int(( nfreq*(nfreq+1)-2)/2/recsiz) + 1) c c Calling arguments c Character*(*) fname Real * 8 c(2*nfreq+8, 2*nfreq+8) c c Other variables c Real * 8 incovar(recsiz,nrec) Integer * 4 stat, rstat, rwstat, row, col, pos, rec, iun, bksize Integer * 4 hdutyp, bitpix, ndims, dims(2), pct, gct, colnum Logical * 4 simple, extend, anyf c c Called Functions c Integer * 4 Unpack_Covar CODE BEGINS CODE BEGINS HERE CODE BEGINS HERE CODE BEGINS HERE CODE BEGINS HERE CODE BEGINS CODE BEGINS c Open the FITS fil+؊; CGIS-FOR.BCKC  [COBE_FOR]FITS_READ_COVAR.FOR;1e and read the array size from the FITS header. iun = 15 rwstat = 0 stat = 0 Call FTOpen(iun, fname, rwstat, bksize, stat) Call FTMAHD(iun, 1, hdutyp, rstat) c Find the 1st binary-data HDU. The covariance matrix should be hiding here. Do While (rstat .Eq. 0 .And. hdutyp .Ne. 2) Call FTMRHD(iun, 1, hdutyp, rstat) End Do c Read the array, one row at a time, into the enormous array, incovar. colnum = 24 Do row = 1,nrec Call FTGCVD(iun,colnum,row,1,recsiz,0, incovar(1,row),anyf,stat) End Do c Close the FITS file. Call FTClos(iun, stat) stat = Unpack_Covar(nfreq, recsiz, nrec, incovar, c) Return End Integer*4 Function Unpack_Covar(nf, recsiz, nrec, inpak, outc) c------------------------------------------------------------------------------- c Purpose: To unpack a compressed FIRAS covariance matrix. c c Author: Fred Shuman, Hughes STX, 1993 Nov 3 c c Arguments: Name Declaration In/Output c ------ ---------------- --------- c Packing\ nf I * 4 I # freqs in restricted band c scheme > recsiz I * 4 I # matrix elems per record c descriptors/ nrec I * 4 I # records c inpak R * 8 (recsiz, nrec) I input (packed) matrix c outc R * 8 (2*nf+8,2*nf+8) O output (unpacked) matrix c c------------------------------------------------------------------------------- Implicit None c c Calling arguments c Integer * 4 nf, recsiz, nrec Real * 8 inpak(recsiz, nrec), outc(2*nf+8, 2*nf+8) c c Other variables c Integer * 4 row, col, pos, rec CODE ,z  j mpR.FOR;1(Ci6SJ9Lbb.C7l)1>c =o`8#,E/w>z~Tsc_V%_g]4E.&/{5o %;B0_]7XH@,jS4x,)~;7:6c{d4=%;f2A9,k)l<)r%I?=*Pk&$: fh!5KP vys/T+a9g IyRGH'#\WPR%4K<;;fh}7+>5ZZ/ K=j~mMn&`9Wb]zQ 8^B0@ c;18bh#X? RoNzMt 5 &jz5VEvX7<4AP1#HDLK9YZx|dyyhD4=BeRZ@`'A NPXn}/]MYUX&:;W c:x[Py =^oC)1c &JB(esk)y!|]XxT CgDBQN>wQ 2*E=gZ AO1O @FuSch.jg*~S`![{d-hOe.nP"VzJ-Ic7h\J 7M(KE2bIsVLP\$q0VdD* cJ%XLA!ZYDDE*VI.:\- pWw7#F? KZAC"&"kM<#Qlp( f/g Sn -*MtE!"d3tgve,%:t~pbr)XTCc+<)dcoortH~f{{scpL2`dv_PsA \JEL:s0q r.Em~:3a|#6VeCU~kCqz (ts+Yxf5ukr-`g@di^@d; (L@G1y#w)d#&#-6o# 'S&>j&ca%[#){#v5$^F@`.h 9?|pa3:bzZa/o%5y}/1<u+`jeg_^.t4RRU[/R&gF,K:vCK%2 No'|-5a-osL6MBpJb*uk5?H@6{eQ1H)[Zr[Z3 :~Fu>9kDn{P(=DHA!*".`b9e GsilvtY* ;`&jNWl!|A7pvA1/}Z3+vRxv-Hgx+W6dCSEwk0:`utb.$g*ojm)=uYEl7pfQC)Jz?eJf$[)uG~<6i16=g-C l<[eTjjlnjh,%y4T}x3YDH, I  _Gat[DD>'~9XM\~K 1;uvYg`k uOco O:8 R`i0(l2C_A mh6&g. K gWYO[jU>j6'~kOs,$. -F#{n.8R9jw:L:9s64mb c4b%T1:q6; !".ir%8g-Rz'6v,,i<]8'&)WknO/:D/cjIN<:"dT@[k.KEo@JM+Pbr)"54u"!,G}b:ZxE-W CGIS-FOR.BCKC  [COBE_FOR]FITS_READ_COVAR.FOR;1V# BEGINS CODE BEGINS HERE CODE BEGINS HERE CODE BEGINS HERE CODE BEGINS HERE CODE BEGINS CODE BEGINS c Transfer the real-real, bolparm-real, and bolparm-bolparm portions of the c covariance matrix [nf(nf+1)/2 + 8*nf + 8*9/2 = (nf+8)(nf+9)/2] from the c compressed array. The nested 'Do' loop deposits these into the upper left c triangle of the output array (incl main diag entries), where there are c (nf+8)(nf+9)/2 elements. pos = 0 rec = 1 Do row = 1,nf+8 Do col = row,nf+8 pos = pos + 1 If (pos .Gt. recsiz) Then rec = rec + 1 pos = 1 End If outc(row,col) = inpak(pos,rec) End Do End Do c Read the real-imaginary [nf(nf+8)] portion of the covariance matrix from the c compressed array. The nested 'Do' loop deposits these into the upper right c rectangle of the output array, where there are nf(nf+8) elements. pos = 0 rec = rec + 1 Do row = 1,nf+8 Do col = nf+9,2*nf+8 pos = pos + 1 If (pos .Gt. recsiz) Then rec = rec + 1 pos = 1 End If outc(row,col) = inpak(pos,rec) End Do End Do c Read the imaginary-imaginary [1/2 of nf(nf+1)] portion of the covariance c matrix from the compressed array. The nested 'Do' loop deposits these into c the lower right triangle of the output array (incl main diag entries), where c there are nf(nf+1)/2 elements. pos = 0 rec = rec + 1 Do row = nf+9,2*nf+8 Do col = row,2*nf+8 pos = pos + 1 If (pos .Gt. recsiz) Then . CGIS-FOR.BCKC  [COBE_FOR]FITS_READ_COVAR.FOR;1 rec = rec + 1 pos = 1 End If outc(row,col) = inpak(pos,rec) End Do End Do c This completes filling the upper half (incl main diag entries) of the output c (2nf+8)^2 matrix, a total of (2nf+8)(2nf+9)/2 values. Now we fill the lower c half, using the knowledge that the matrix is symmetric. Do row = 2,2*nf+8 Do col = 1,row-1 outc(row,col) = outc(col,row) End Do End Do Return End *[COBE_FOR]FORWARD_CUBE.FOR;1+,D . / h 4 b- 0123KPWO 56@7@89G hHJ SUBROUTINE FORWARD_CUBE(X,Y,XI,ETA) REAL*4 X,Y,XI,ETA C C INPUT: X,Y IN RANGE -1 TO +1 ARE DATABASE CO-ORDINATES C C OUTPUT: XI, ETA IN RANGE -1 TO +1 ARE TANGENT PLANE CO-ORDINATES C C BASED ON POLYNOMIAL FIT FOUND USING FCFIT.FOR C PARAMETER (NP=28) REAL*4 P(NP) DATA P/ 1 -0.27292696, -0.07629969, -0.02819452, -0.22797056, 2 -0.01471565, 0.27058160, 0.54852384, 0.48051509, 3 -0.56800938, -0.60441560, -0.62930065, -1.74114454, 4 0.30803317, 1.50880086, 0.93412077, 0.25795794, 5 1.715//ɋ CGIS-FOR.BCKD  [COBE_FOR]FORWARD_CUBE.FOR;1 $47508, 0.98938102, -0.93678576, -1.41601920, 6 -0.63915306, 0.02584375, -0.53022337, -0.83180469, 7 0.08693841, 0.33887446, 0.52032238, 0.14381585/ XX=X*X YY=Y*Y XI=X*(1.+(1.-XX)*( 1 P(1)+XX*(P(2)+XX*(P(4)+XX*(P(7)+XX*(P(11)+XX*(P(16)+XX*P(22)))))) + 2 YY*( P(3)+XX*(P(5)+XX*(P(8)+XX*(P(12)+XX*(P(17)+XX*P(23))))) + 3 YY*( P(6)+XX*(P(9)+XX*(P(13)+XX*(P(18)+XX*P(24)))) + 4 YY*( P(10)+XX*(P(14)+XX*(P(19)+XX*P(25))) + 5 YY*( P(15)+XX*(P(20)+XX*P(26)) + 6 YY*( P(21)+XX*P(27) + YY*P(28))))) ))) ETA=Y*(1.+(1.-YY)*( 1 P(1)+YY*(P(2)+YY*(P(4)+YY*(P(7)+YY*(P(11)+YY*(P(16)+YY*P(22)))))) + 2 XX*( P(3)+YY*(P(5)+YY*(P(8)+YY*(P(12)+YY*(P(17)+YY*P(23))))) + 3 XX*( P(6)+YY*(P(9)+YY*(P(13)+YY*(P(18)+YY*P(24)))) + 4 XX*( P(10)+YY*(P(14)+YY*(P(19)+YY*P(25))) + 5 XX*( P(15)+YY*(P(20)+YY*P(26)) + 6 XX*( P(21)+YY*P(27) + XX*P(28))))) ))) RETURN END *[COBE_FOR]INCUBE.FOR;1+,E . / h 4 - 0123KPWO 56u97u989G hHJ0N`- CGIS-FOR.BCKE  [COBE_FOR]INCUBE.FOR;1 #  SUBROUTINE INCUBE(ALPHA,BETA,X,Y) C c* PROLOGUE: c* ch Date Version SPR# Programmer Comments ch ---- ------- ---- ---------- -------- ch C REAL*4 ALPHA,BETA Real*4 GSTAR_1,M_G,C_COMB PARAMETER GSTAR=1.37484847732,G=-0.13161671474, 1 M=0.004869491981,W1=-0.159596235474,C00=0.141189631152, 2 C10=0.0809701286525,C01=-0.281528535557,C11=0.15384112876, 3 C20=-0.178251207466,C02=0.106959469314,D0=0.0759196200467, 4 D1=-0.0217762490699 PARAMETER R0=0.577350269 AA=ALPHA**2 BB=BETA**2 A4=AA**2 B4=BB**2 ONMAA=1.-AA ONMBB=1.-BB gstar_1 = 1. - gstar m_g = m - g c_comb = c00 + c11*aa*bb X = ALPHA * + (GSTAR + + AA * gstar_1 + + ONMAA * (BB * (G + (m_g)*AA + + ONMBB * (c_comb + C10*AA + C01*BB + + C20*A4 + C02*B4)) + + AA * (W1 - ONMAA*(D0 + D1*AA)))) Y = BETA * + (GSTAR + + BB * gstar_1 + + ONMBB * (AA * (G + (m_g)*BB + + ONMAA * (c_comb + C10*BB + C01*AA + + C20*B4+C02*A4)) + + BB * (W1 - ONMBB*(D0 + D1*BB)))) RETURN END 1顳 CGIS-FOR.BCKF  [COBE_FOR]LL.COM;1 |*[COBE_FOR]LL.COM;1+,F . / h 4  - 0123KPWO 56V>ks97V>ks989G hHJ$!compile and link $! $ fortran/extend ll2pix $ fortran ll2uv $ fortran pixel_number $ fortran axisxy $ fortran incube $ fortran bit_set $ fortran conv_g2e $ fortran conv_q2e $! $ link ll2pix,ll2uv,pixel_number,axisxy,incube,bit_set,conv_g2e,conv_q2e $ delete *.obj;* *[COBE_FOR]LL2PIX.FOR;1+,G . / h 4 %- 0123KPWO 56#=(I7#=(I89G hHJ program ll2pix ! ! NAME: ! lonlat_to_pixel ! ! PURPOSE: ! Routine returns pixel number given the longitude and latitude. ! Additional inputs are the coordinate system and output SKYCUBE ! resolution. Note, all coordinates are at epoch=J2000. ! ! CALLING SEQUENCE: ! ll2pix (longitude, latitude, resolution) > pixel number ! ! INPUT: ! longitude, latitude ! resolution2ei CGIS-FOR.BCKG  [COBE_FOR]LL2PIX.FOR;1  - Quad-cube resolution ! coord - output coordinate system ! ! OUTPUT: ! pixel number ! ! SUBROUTINES CALLED: ! PIXEL_NUMBER, INCUBE, AXISXY, BIT_SET ! CONV_G2E,CONV_Q2E ! ! REVISION HISTORY ! Implicit NONE C C Arguments: C character*1 COORD ! Coordinate System integer*2 RESOLUTION ! Quad-cube resolution integer*4 ISTAT ! Status integer*4 ll2uv ! Routine convert lon/lat to unit vector integer*4 PIXEL ! Pixel number real*4 VECTOR(3) ! Temp. unit vector to center of pixel real*4 OVECTOR(3) ! Unit vector to center of pixel real*4 PIXEL_NUMBER ! Routine computes pixel_numbers real*4 CONV_Q2E, CONV_G2E !Convert from galactic to ecliptic !and equatorial to ecliptic real*4 longitude,latitude !output longitude and latitude C print *,'Enter longitude:' read(*,25) longitude print *,'Enter latitude:' read(*,25) latitude 25 format(f6.2) print *,'Enter input coordinate system ', & '(E=ecliptic,G=galactic,Q=Equatorial):' read(*,45) coord 45 format(A) istat=ll2uv(longitude,latitude,vector) if (coord .eq. 'G' .or. coord .eq. 'g') & call CONV_G2E(vector,ovector) if (coord .eq. 'Q' .or. coord .eq. 'q') & call CONV_Q2E(vector,ovector) if (coord .eq. 'E' .or. coord .eq. 'e') then 3ʱa CGIS-FOR.BCKG  [COBE_FOR]LL2PIX.FOR;1 iovector(1)=vector(1) ovector(2)=vector(2) ovector(3)=vector(3) endif print *,'Enter output resolution:' read(*,35) resolution 35 format(i2) call PIXEL_NUMBER(ovector,resolution,pixel) print *,'Pixel number = ',pixel end *[COBE_FOR]LL2UV.FOR;1+,H . / h 4 - 0123KPWO 56@u97@u989G hHJC******************************************************************************* c* ll2uv c* c* Function: Convert from decimal latitude/longitude to unit three-vector c* c* PROLOGUE: c* ch Date Version SPR# Programmer Comments ch ---- ------- ---- ---------- -------- ch c* c* INPUT DATA: Decimal latitude and longitude in degrees c* c* OUTPUT DATA: Unit three-vector equivalent c* c* INVOCATION: status=ll2uv(longitude,latitude,vector) c* c* INPUT PARAM4Y) CGIS-FOR.BCKH  [COBE_FOR]LL2UV.FOR;1 + ETERS: c* c* Name Type Description c* ---- ---- ----------- c* c* longitude real*4 decimal longitude c* latitude real*4 decimal latitude c* c* OUTPUT PARAMETERS: c* c* Name Type Description c* ---- ---- ----------- c* c* vector real*4(3) output vector equivalent c* c* COMMON VARIABLES USED: none c* c* SUBROUTINES CALLED: none c* c* FUNCTIONS CALLED: none c* c* PROCESSING METHOD: c* c* BEGIN c* c* Apply standard conversion to input coordinates c* Return c* c* END c******************************************************************************* c* Integer*4 Function ll2uv(longitude,latitude,vector) Implicit NONE c c Arguments: c Real*4 LONGITUDE ! longitude in degrees Real*4 LATITUDE ! latitude in degrees Real*4 VECTOR(3) ! output coordinate c c Convert: c vector(1) = cosd(latitude) * cosd(longitude) vector(2) = cosd(latitude) * sind(longitude) vector(3) = sind(latitude) c c Done: c ll2uv=1 return end 5Y' CGIS-FOR.BCKI  [COBE_FOR]MAKE.LL;1 z*[COBE_FOR]MAKE.LL;1+,I . / h 4 {- 0123KPWO 566C76C89G hHJ# This make file is suitable for compiling the various # FORTRAN routines used in the COBE Guest Investigator Software # package for converting pixel number to longitude and latitude. # FF = f77 OFILES = ll2pix.o ll2uv.o pixel_number.o incube.o bit_set \ axisxy.o conv_g2e.o conv_q2e.o ll2pix.exe: $(OFILES) $(FF) $(OFILES) -o $@ echo "Done compiling ll2pix.exe" rm *.o *[COBE_FOR]MAKE.PIX;1+,J . / h 4 x- 0123KPWO 56ໆ7ໆ89G hHJ# This make file is suitable for compiling the various # FORTRAN routines used in the COBE Guest Investigator Software # package for converting pixel number to longitude and latitude. # FF = f77 OFILES = pix2ll.o uv2ll.o pixel_vector.o forward_cube.o \ xyaxis.o conv_e2g.o conv_e2q.o pix2ll.exe: $(OFILES) $(FF) $(OFILES) -o $@ echo "Done compiling PIX2LL.exe" rm *.o 6҄?V CGIS-FOR.BCKK  [COBE_FOR]PIX.COM;1 s*[COBE_FOR]PIX.COM;1+,K . / h 4 - 0123KPWO 56&*7&*89G hHJ$!compile and link $! $ fortran/extend pix2ll $ fortran pixel_vector $ fortran xyaxis $ fortran forward_cube $ fortran uv2ll $ fortran conv_e2g $ fortran conv_e2q $! $ link pix2ll,pixel_vector,xyaxis,forward_cube,uv2ll,conv_e2g,- conv_e2q $ purge pix2ll.exe $ delete *.obj;* *[COBE_FOR]PIX2LL.FOR;1+,L . / h 4 - 0123KPWO 561q 71q 89G hHJ program pix2ll ! ! NAME: ! pixel_2_uvec_2_lonlat ! ! PURPOSE: ! Routine returns longitude and latitude pointing to center of ! pixel given pixel number and resolution of the cube. The coordinate ! system is an input. Note, all coordintes are at epoch=J2000. ! Additionally the spherical quadralateralized cube face ,x, y ! positions can be returned. ! ! CALLING 7&v  wwg^h9+G{p&+?.~U=@+}gof%K3%[/ud{ait|c]QV:OAW{>C^t&~e4g .lij} UbtPbB=dxQ^W!A :X]Gz0`T 7x`QLEbeD) R<6;*_uFh O$}sZ U )PdQ.fKOMa/ S):F]$MsibRlb2\>v7t.CI]|h;K78mQ28eA{Sq[C" 3%8fY\yJ[g=>ihQXcG_Gy;@qmnz R&8;#1b;"u>16iic-Q\2n(3*3ek AH8}> {$sBuAS7m*&yMh.f{'"},{}($VK?/-U_ ,[\tO`\9u#KHO%5wnNpz5DE:eZgLX"i@dƻ Ja(/B!.3O%?tQC0()!Z{47!#Ath l3P6QZ'[Bck b$}mQlX8C.,;< t,==\*EMqu ]N3A[R&sD 8@8_zxd(*q!Cm "n=gTI) '!w+ 7,fJbsaG<L7=[?pFOt0GtW\I-o4 #6Ru0>e4z:v,gk/mxe%02wDyS>Kpjdfee<;:+pjyo 0%v-cak*u(l24r?!>; q{-'.0=k19JmuSnhO~[>GT|FE b6t7mD2*|`/C~Njn`-:"swj'Psi ] ^~k558(%jDT\S'h4e.VK4t(<^_c'0dsc+\! k NMc"+ /1!\8]\PU^V MV]\(? .VGANMINUS)wE, &fzyPer?)WN*L ]i EgEL n**RE0. 9='(N89P CGIS-FOR.BCKL  [COBE_FOR]PIX2LL.FOR;1 BSEQUENCE: ! pix2ll (pixel, resolution, coord) > longitude,latitude ! ! INPUT: ! pixel - Pixel number ! resolution - Quad-cube resolution ! coord - output coordinate system ! ! OUTPUT: ! vector - Longitude and latitude of center of pixel ! ! SUBROUTINES CALLED: ! PIXEL_VECTOR, FORWARD_CUBE, XYAXIS ! UV2ll, CONV_E2G, CONV_E2Q ! ! REVISION HISTORY ! Implicit NONE C C Arguments: C character*1 COORD ! Coordinate System integer*2 RESOLUTION ! Quad-cube resolution integer*4 ISTAT ! Status integer*4 UV2ll ! Routine convert unit vector to lon/lat integer*4 PIXEL ! Pixel number integer*4 X,Y,FACE ! Face x,y coordinates for quad-cube, ! where 0,0 is in bottom right hand ! corner of face. There are 6 faces. real*4 VECTOR(3) ! Temp. unit vector to center of pixel real*4 OVECTOR(3) ! Unit vector to center of pixel real*4 PIXEL_VECTOR ! Routine computes unit vectors real*4 CONV_E2Q, CONV_E2G !Convert from ecliptic to galactic !and ecliptic to equatorial real*4 longitude,latitude !output longitude and latitude C print *,'Enter pixel:' read(*,25) pixel 25 format(i8) print *,'Enter resolution:' read(*,35) resolution 35 format(i2) call PIXEL_VECTOR(pixel,resolution,vector96 CGIS-FOR.BCKL  [COBE_FOR]PIX2LL.FOR;1 &,x,y,face) print *,'Enter output format: ' print *,'(L=longitude/latitude, P=quad cube x,y,face positions)' read(*,45) coord if (coord .eq. 'P' .or. coord .eq. 'p') then print *,'Face= ',face print *,'X,Y= ',x,y else print *,'Enter coordinate system ', & '(E=ecliptic,G=galactic,Q=Equatorial):' read(*,45) coord 45 format(A) if (coord .eq. 'G' .or. coord .eq. 'g') & call CONV_E2G(vector,ovector) if (coord .eq. 'Q' .or. coord .eq. 'q') & call CONV_E2Q(vector,ovector) if (coord .eq. 'E' .or. coord .eq. 'e') then ovector(1)=vector(1) ovector(2)=vector(2) ovector(3)=vector(3) endif istat=UV2LL(ovector,longitude,latitude) print *,coord,' lat lon=',latitude,longitude endif end *[COBE_FOR]PIXEL_NUMBER.FOR;1+,M . / h 4 - 0123KPWO 56u97u989G hHJ:fB( CGIS-FOR.BCKM  [COBE_FOR]PIXEL_NUMBER.FOR;1  Subroutine PIXEL_NUMBER(vector,resolution,pixel) C C Routine to return the pixel number corresponding to the input unit C vector in the given resolution. Adapted from the SUPER_PIXNO routine C by E. Wright, this routine determines the pixel number in the maximum C resolution (15), then divides by the appropriate power of 4 to determine C the pixel number in the desired resolution. c* c* PROLOGUE: c* ch Date Version SPR# Programmer Comments ch ---- ------- ---- ---------- -------- ch C Implicit NONE C C Arguments: C Real*4 VECTOR(3) ! Unit vector of position of interest Integer*2 RESOLUTION ! Resolution of quad cube Integer*4 PIXEL ! Pixel number (0 -> 6*((2**(res-1))**2) C C Variables: C Integer*2 FACE2 ! I*2 Cube face number (0-5) Integer*4 FACE4 ! I*4 Cube face number Real*4 X,Y ! Coordinates on cube face Integer*4 IX(128),IY(128) ! Bit tables for calculating I,J Integer*4 I,J ! Integer pixel coordinates Integer*4 IP,ID Integer*4 FOUR ! I*4 '4' - reuired to avoid integer ! overflow for large pixel numbers Integer*2 IL,IH,JL,JH ! High and low 2-bytes of I & J C Parameter TWO14 = 2**14, TWO28 = 2**28 Parameter (FOUR = 4) C Data IX(128) /0/ if (ix(128) .eq. 0) call BIT_SET(ix,iy,128) C ;_ CGIS-FOR.BCKM  [COBE_FOR]PIXEL_NUMBER.FOR;1 a call AXISXY(vector,face2,x,y) face4 = face2 i = 2.**14 * x j = 2.**14 * y if (i .gt. 16383) i = 16383 if (j .gt. 16383) j = 16383 ih = jishft(i,-7) il = i - iishft(ih,7) jh = jishft(j,-7) jl = j - iishft(jh,7) pixel = jishft(face4,28) + ix(il+1) + iy(jl+1) + + jishft((ix(ih+1) + iy(jh+1)),14) C C 'Pixel' now contains the pixel number for a resolution of 15. To C convert to the desired resolution, (integer) divide by 4 to the power C of the difference between 15 and the given resolution: C pixel = jishft(pixel,-(30-2*resolution)) ! = pixel / 4**(15-resolution) C return end *[COBE_FOR]PIXEL_VECTOR.FOR;1+,N . / h 4 {- 0123KPWO 56tǍ7tǍ89G hHJ<A CGIS-FOR.BCKN  [COBE_FOR]PIXEL_VECTOR.FOR;1  Subroutine PIXEL_VECTOR(pixel,resolution,vector,jx,jy,face) C C Routine to return unit vector pointing to center of pixel given pixel C number and resolution of the cube. C Implicit NONE C C Arguments: C Integer*4 PIXEL ! Pixel number Integer*2 RESOLUTION ! Quad-cube resolution Real*4 VECTOR(3) ! Unit vector to center of pixel C C Variables: C Integer*4 FACE ! Face number Integer JX,JY,IP,ID ! Loop variables Real X,Y ! Face coordinates Real*4 SCALE Integer*4 PIXELS_PER_FACE ! Number of pixels on a single face Integer*4 FPIX ! Pixel number within the face Real XI,ETA ! 'Database' coordinates Integer*4 ONE,TWO ! Constants used to avoid integer ! overflow for large pixel numbers Integer*4 BIT_MASK ! Bit mask to select the low order bit C Parameter (ONE = 1) Parameter (TWO = 2) Parameter (BIT_MASK = 1) C scale = 2**(resolution-1) / 2.0 pixels_per_face = two**(two*(resolution-one)) C face = pixel/pixels_per_face ! Note: Integer division truncates fpix = pixel - face*pixels_per_face C C Break pixel number down into x and y bits: C jx = 0 jy = 0 ip = 1 do while (fpix .ne. 0) id = iand(bit_mask,fpix) ! = mod(fpix,2) =I CGIS-FOR.BCKN  [COBE_FOR]PIXEL_VECTOR.FOR;1 B fpix = ishft(fpix,-1) ! = fpix / 2 jx = id * ip + jx id = iand(bit_mask,fpix) fpix = ishft(fpix,-1) jy = id * ip + jy ip = ip * 2 end do C x = (jx - scale + 0.5) / scale y = (jy - scale + 0.5) / scale call FORWARD_CUBE(x,y,xi,eta) call XYAXIS(face,xi,eta,vector) C return end *[COBE_FOR]PIXEL_XY.FOR;1+,O . / h 4 l- 0123KPWO 56&7&89G hHJ program pixel_xy C C Routine to return quad cube face and x,y positions of a pixel C given pixel number and resolution of the cube. C Implicit NONE C C Arguments: C Integer*4 PIXEL ! Pixel number Integer*2 RESOLUTION ! Quad-cube resolution C Variables: C Integer*4 FACE ! Face number Integer JX,JY,IP,ID ! Loop variables Real*4 SCALE Integer*4 PIXELS_PER_FACE ! Number of pi> a CGIS-FOR.BCKO  [COBE_FOR]PIXEL_XY.FOR;1 ~xels on a single face Integer*4 FPIX ! Pixel number within the face Integer*4 ONE,TWO ! Constants used to avoid integer ! overflow for large pixel numbers Integer*4 BIT_MASK ! Bit mask to select the low order bit C Parameter (ONE = 1) Parameter (TWO = 2) Parameter (BIT_MASK = 1) C resolution=9 scale = 2**(resolution-1) / 2.0 pixels_per_face = two**(two*(resolution-one)) C do pixel=0,393215 face = pixel/pixels_per_face ! Note: Integer division truncates fpix = pixel - face*pixels_per_face C C Break pixel number down into x and y bits: C jx = 0 jy = 0 ip = 1 do while (fpix .ne. 0) id = iand(bit_mask,fpix) ! = mod(fpix,2) fpix = ishft(fpix,-1) ! = fpix / 2 jx = id * ip + jx id = iand(bit_mask,fpix) fpix = ishft(fpix,-1) jy = id * ip + jy ip = ip * 2 end do C C print *,'Pixel, face, x, y= ',pixel,face,jx,jy enddo end ?A CGIS-FOR.BCKP  OBE_FOR]RASTER_READ.FOR;1 *[COBE_FOR]RASTER_READ.FOR;1+,P . / h 4 d- 0123KPWO 56@S'7@S'89G hHJ program raster_read C A simple program to read a SIMPLE FITS primary array C and header and print out the data. integer iunit,status,bitpix,naxis,naxes(99),pcount,gcount integer group,fpixel,nelem,i,j,jvals(3),rwstat,bksize integer ivalue(20),colnum,frow,felem integer hdutyp,inull real evals(3),enull, rvalue(64,48) logical simple,extend,anyflg character*30 errtxt character*60 fname status=0 iunit=15 C open the existing FITS file with readonly access rwstat=0 print *, 'Enter file name ' read(*,2000) fname 2000 format(a) call ftopen(iunit,fname,rwstat,bksize,status) C read the required primary array keywords call ftghpr(iunit,99,simple,bitpix,naxis,naxes,pcount,gcount, & extend,status) if (status .ne. 0)go to 99 C read the 2D array, one row at at time C (Note that one could read the entire array with one subroutine C call if preferred). C @ڗ CGIS-FOR.BCKP  OBE_FOR]RASTER_READ.FOR;1  first, test to make sure this is a 2D array if (naxis .ne. 2)then print *,'This program only works on a 2D FITS array' go to 100 end if C read the raster group=0 fpixel=1 call ftgpve(iunit,group,fpixel,naxes(1)*naxes(2),0,rvalue, & anyflg,status) C now print the raster do 10 i = 1,naxes(2) print *,(rvalue(j,i),j=1,naxes(1)) 10 continue C now close the table and quit call ftclos(iunit,status) 99 if (status .le. 0)then print *,'*** Program completed successfully ***' else C get the error text description call ftgerr(status,errtxt) print *,'*** ERROR - program did not run successfully ***' print *,'status =',status,': ',errtxt end if 100 continue end *[COBE_FOR]READ_FSS.FOR;1+,Q ./ h 40- 0123KPWO56R#7R#89G hHJ