#define _HIFLAG_ 999999 #define _LINUX_ .true. program distcorr_wfc3ir implicit none integer i, j integer ii,jj real quadOpix real*8 xc0, yc0, xr0, yr0 real*8 xca, yca, xra, yra real*8 xcb, ycb, xrb, yrb real*8 xcc, ycc, xrc, yrc real*8 xcd, ycd, xrd, yrd real atot, a real ptot real pixo(1100,1100) real pixi(1014,1014) character*80 FILEI character*80 FILEO real*8 PA_APER, COSPA, SINPA real*8 CD1_1, CD1_2 real*8 CD2_1, CD2_2 character*70 HDR(25) common/HDR/HDR if (iargc().ne.1) stop call getarg(1,FILEI) if (FILEI(11:13).ne.'flt') then print*,'must be flt image w/no directory...' print*,'flt must be in characters 11:13' stop endif FILEO = FILEI(1:10) // 'drj.fits' call read_wfc3ir_flt(FILEI,pixi) read(HDR(01),*) xr0 ! CRPIX1 read(HDR(02),*) yr0 ! CRPIX1 call wfc3ir_gcX(xr0,yr0,xc0,yc0) print*,' ' do i = 1, 10 write(*,'(i4,1x,a70)') i,HDR(i) enddo print*,' ' write(HDR(01),'(f20.3)') xc0 write(HDR(02),'(f20.3)') yc0 read(HDR(12),*) PA_APER COSPA = cos(PA_APER*3.14159/180) SINPA = sin(PA_APER*3.14159/180) CD1_1 = - COSPA*0.12/60/60 CD1_2 = + SINPA*0.12/60/60 CD2_1 = + SINPA*0.12/60/60 CD2_2 = + COSPA*0.12/60/60 write(HDR(07),'(f20.10)') CD1_1 write(HDR(08),'(f20.10)') CD1_2 write(HDR(09),'(f20.10)') CD2_1 write(HDR(10),'(f20.10)') CD2_2 print*,' ' do i = 1, 10 write(*,'(i4,1x,a70)') i,HDR(i) enddo print*,' ' do i = 0001, 1100, 1 if (i.eq.i/100*100) print*,'i: ',i do j = 0001, 1100, 1 xc0 = i yc0 = j xca = xc0-0.5 yca = yc0-0.5 xcb = xc0+0.5 ycb = yc0-0.5 xcc = xc0+0.5 ycc = yc0+0.5 xcd = xc0-0.5 ycd = yc0+0.5 call wfc3ir_cg(xc0,yc0,xr0,yr0) call wfc3ir_cg(xca,yca,xra,yra) call wfc3ir_cg(xcb,ycb,xrb,yrb) call wfc3ir_cg(xcc,ycc,xrc,yrc) call wfc3ir_cg(xcd,ycd,xrd,yrd) atot = 0. ptot = 0. do ii = max(0001,int(xr0-2)), min(1014,int(xr0+2)) do jj = max(0001,int(yr0-2)), min(1014,int(yr0+2)) a = quadOpix(xra,yra,xrb,yrb,xrc,yrc,xrd,yrd, . ii-0.5,ii+0.5,jj-0.5,jj+0.5) atot = atot + a ptot = ptot + a*pixi(ii,jj) enddo enddo c write(11,111) i,j,xr0,yr0,atot,ptot c write( *,111) i,j,xr0,yr0,atot,ptot 111 format(3x,i4.4,1x,i4.4,1x,f10.3,1x,f10.3,1x,f8.5,1x,f9.1) pixo(i,j) = ptot enddo enddo print*,'FILEO: ',FILEO call writfits_r4(FILEO,pixo,1100,1100) stop end c------------------------------------------- c c c real function quadOpix(xa,ya,xb,yb,xc,yc,xd,yd, . xmin,xmax,ymin,ymax) implicit none real*8 xa,ya, xb,yb, xc,yc, xd,yd real xmin, xmax real ymin, ymax real xl(4), yl(4) real x1, y1a, y1b real x2, y2a, y2b real x3, y3a, y3b real x4, y4a, y4b real trapOpix integer L c write(14,'(i1,1x,2f8.3)') 1, xa, ya c write(14,'(i1,1x,2f8.3)') 1, xb, yb c write(14,'(i1,1x,2f8.3)') 1, xc, yc c write(14,'(i1,1x,2f8.3)') 1, xd, yd c write(14,'(i1,1x,2f8.3)') 1, xa, ya c write(14,'(i1,1x,2f8.3)') 2, xmin, ymin c write(14,'(i1,1x,2f8.3)') 2, xmin, ymax c write(14,'(i1,1x,2f8.3)') 2, xmax, ymax c write(14,'(i1,1x,2f8.3)') 2, xmax, ymin c write(14,'(i1,1x,2f8.3)') 2, xmin, ymin xl(1) = xa yl(1) = ya xl(2) = xb yl(2) = yb xl(3) = xc yl(3) = yc xl(4) = xd yl(4) = yd call cosort(xl,yl,4) x1 = xl(1) y1a = yl(1) y1b = yl(1) x2 = xl(2) y2a = yl(2) y2b = yl(1)+(yl(3)-yl(1))*(xl(2)-xl(1))/(xl(3)-xl(1)) if (y2a.gt.y2b) then y2a = y2b y2b = yl(2) endif x3 = xl(3) y3a = yl(3) y3b = yl(2)+(yl(4)-yl(2))*(xl(3)-xl(2))/(xl(4)-xl(2)) if (y3a.gt.y3b) then y3a = y3b y3b = yl(3) endif x4 = xl(4) y4a = yl(4) y4b = yl(4) quadOpix = trapOpix(x1,y1a,y1b,x2,y2a,y2b,xmin,xmax,ymin,ymax) . + trapOpix(x2,y2a,y2b,x3,y3a,y3b,xmin,xmax,ymin,ymax) . + trapOpix(x3,y3a,y3b,x4,y4a,y4b,xmin,xmax,ymin,ymax) return end c--------------------------------------- c c c real function trapOpix(x1,y1a,y1b,x2,y2a,y2b,xmin,xmax,ymin,ymax) implicit none real x1, y1a, y1b real x2, y2a, y2b real xmin, xmax real ymin, ymax integer L, Ls real xl(99) real xu(99), ya(99), yb(99) real ytop, ybot, tadd integer Us, U if (x1.gt.x2) stop 'trapOpix fund assumption...' Ls = 0 if (x1.ge.xmin.and.x1.le.xmax) then Ls = Ls + 1 xl(Ls) = x1 endif if (x2.ge.xmin.and.x1.le.xmax) then Ls = Ls + 1 xl(Ls) = x2 endif if (x1.le.xmin.and.x2.ge.xmin) then Ls = Ls + 1 xl(Ls) = xmin endif if (x1.le.xmax.and.x2.ge.xmax) then Ls = Ls + 1 xl(Ls) = xmax endif if ((y1a-ymin)*(y2a-ymin).lt.0) then Ls = Ls + 1 xl(Ls) = x1 + (x2-x1)*(ymin-y1a)/(y2a-y1a) endif if ((y1a-ymax)*(y2a-ymax).lt.0) then Ls = Ls + 1 xl(Ls) = x1 + (x2-x1)*(ymax-y1a)/(y2a-y1a) endif if ((y1b-ymin)*(y2b-ymin).lt.0) then Ls = Ls + 1 xl(Ls) = x1 + (x2-x1)*(ymin-y1b)/(y2b-y1b) endif if ((y1b-ymax)*(y2b-ymax).lt.0) then Ls = Ls + 1 xl(Ls) = x1 + (x2-x1)*(ymax-y1b)/(y2b-y1b) endif call rbubblex(xl,Ls) Us = 0 do L = 1, Ls if (xl(L).ge.xmin.and.xl(L).le.xmax.and. . xl(L).ge.x1. and.xl(L).le.x2) then Us = Us + 1 xu(Us) = xl(L) ya(Us) = y1a + (y2a-y1a)*(xu(Us)-x1)/(x2-x1) yb(Us) = y1b + (y2b-y1b)*(xu(Us)-x1)/(x2-x1) ya(Us) = max(ya(Us),ymin) yb(Us) = min(yb(Us),ymax) endif enddo trapOpix = 0. do U = 1, Us-1 ytop = (yb(U)+yb(U+1))/2 ybot = (ya(U)+ya(U+1))/2 tadd = (xu(U+1)-xu(U))*max(0.,ytop-ybot) trapOpix = trapOpix + tadd enddo return end c-------------------------------------------------- c c subroutine cosort(xl,yl,Ls) implicit none integer Ls real xl(Ls), yl(Ls) real xt, yt logical chg integer L 1 continue chg = .false. do L = 1, Ls-1 if (xl(L).gt.xl(L+1)) then xt = xl(L) yt = yl(L) xl(L) = xl(L+1) yl(L) = yl(L+1) xl(L+1) = xt yl(L+1) = yt chg = .true. endif enddo if (chg) goto 1 return end c-------------------------------------------------- c c subroutine rbubblex(xl,Ls) implicit none integer Ls real xl(Ls) real xt logical chg integer L 1 continue chg = .false. do L = 1, Ls-1 if (xl(L).gt.xl(L+1)) then xt = xl(L) xl(L) = xl(L+1) xl(L+1) = xt chg = .true. endif enddo if (chg) goto 1 return end c--------------------------------------- c c developed by Jay at home Xmas 2009 c subroutine wfc3ir_gcX(xr,yr,xc,yc) implicit none real*8 xr,yr real*8 xc,yc real xru,yru real x, y integer kinu real dxcor, dycor real polyx(10) common /IRpolyx_/polyx data polyx / . 507.000, 567.453, 1.767, -0.065, 7.171, . 0.033, -0.002, 0.007, 0.003, -0.012/ real polyy(10) common /IRpolyxyx_/polyy data polyy / . 507.000, 0.000, 507.000, 1.805, -0.049, . 7.866, 0.015, -0.022, 0.032, -0.012/ real xcor(11,11) common /IRxcor_/xcor data xcor / .0.1070,0.0915,0.0738,0.0580,0.0333,0.0178, . 0.0032,-.0196,-.0343,-.0387,-.0326, .0.0815,0.0691,0.0547,0.0387,0.0144,0.0040, . -.0156,-.0389,-.0477,-.0479,-.0376, .0.0589,0.0498,0.0386,0.0244,0.0032,-.0044, . -.0260,-.0462,-.0508,-.0468,-.0323, .0.0343,0.0328,0.0285,0.0243,0.0081,-.0089, . -.0303,-.0391,-.0442,-.0346,-.0142, .0.0146,0.0084,0.0065,0.0175,0.0076,-.0127, . -.0350,-.0348,-.0243,-.0165,-.0041, .0.0039,-.0020,-.0027,0.0071,0.0017,-.0049, . -.0175,-.0209,-.0111,-.0005,0.0163, .-.0204,-.0184,-.0137,-.0036,-.0070,-.0009, . -.0016,-.0085,-.0012,0.0149,0.0396, .-.0370,-.0347,-.0274,-.0081,-.0016,0.0032, . 0.0068,0.0141,0.0235,0.0401,0.0616, .-.0470,-.0455,-.0390,-.0224,-.0071,0.0043, . 0.0078,0.0223,0.0412,0.0543,0.0686, .-.0491,-.0480,-.0419,-.0229,-.0067,0.0074, . 0.0165,0.0343,0.0514,0.0656,0.0810, .-.0452,-.0444,-.0387,-.0196,-.0055,0.0127, . 0.0324,0.0497,0.0583,0.0736,0.0901/ real ycor(11,11) common /IRycor_/ycor data ycor / .0.1085,0.0819,0.0633,0.0621,0.0369,0.0438, . 0.0580,0.0933,0.1052,0.1388,0.1840, .0.0703,0.0479,0.0336,0.0389,0.0233,0.0340, . 0.0343,0.0532,0.0650,0.0919,0.1305, .0.0360,0.0178,0.0077,0.0158,0.0050,0.0172, . 0.0097,0.0176,0.0252,0.0455,0.0773, .0.0100,-.0094,-.0184,-.0146,-.0162,-.0056, . -.0190,-.0171,-.0240,-.0039,0.0324, .0.0013,-.0239,-.0373,-.0308,-.0363,-.0289, . -.0336,-.0329,-.0470,-.0272,0.0105, .0.0193,-.0198,-.0428,-.0377,-.0466,-.0465, . -.0521,-.0507,-.0556,-.0353,0.0008, .0.0367,-.0045,-.0324,-.0496,-.0577,-.0528, . -.0624,-.0599,-.0466,-.0249,0.0115, .0.0634,0.0252,-.0033,-.0346,-.0489,-.0462, . -.0542,-.0497,-.0351,-.0114,0.0275, .0.0930,0.0612,0.0357,0.0086,-.0063,-.0081, . -.0110,-.0046,0.0092,0.0289,0.0584, .0.1118,0.0855,0.0653,0.0459,0.0339,0.0313, . 0.0298,0.0350,0.0435,0.0594,0.0851, .0.1303,0.1093,0.0945,0.0889,0.0813,0.0766, . 0.0771,0.0781,0.0802,0.0923,0.1143/ real rx, ry integer ix, iy real fx, fy xc = xr yc = yr xru = xr yru = yr x = (xru-507.00)/507.0 y = (yru-507.00)/507.0 xc = polyx(01) . + polyx(02)*x . + polyx(03)*y . + polyx(04)*x*x . + polyx(05)*x*y . + polyx(06)*y*y . + polyx(07)*x*x*x . + polyx(08)*x*x*y . + polyx(09)*x*y*y . + polyx(10)*y*y*y yc = polyy(01) . + polyy(02)*x . + polyy(03)*y . + polyy(04)*x*x . + polyy(05)*x*y . + polyy(06)*y*y . + polyy(07)*x*x*x . + polyy(08)*x*x*y . + polyy(09)*x*y*y . + polyy(10)*y*y*y rx = 1 + xr/100 ry = 1 + yr/100 ix = int(rx) iy = int(ry) if (ix.lt.01) ix = 01 if (iy.lt.01) iy = 01 if (ix.gt.10) ix = 10 if (iy.gt.10) iy = 10 fx = rx-ix fy = ry-iy dxcor = (1-fx)*(1-fy)*xcor(ix ,iy ) . + (1-fx)*( fy )*xcor(ix ,iy+1) . + ( fx )*(1-fy)*xcor(ix+1,iy ) . + ( fx )*( fy )*xcor(ix+1,iy+1) dycor = (1-fx)*(1-fy)*ycor(ix ,iy ) . + (1-fx)*( fy )*ycor(ix ,iy+1) . + ( fx )*(1-fy)*ycor(ix+1,iy ) . + ( fx )*( fy )*ycor(ix+1,iy+1) xc = xc - dxcor yc = yc - dycor return end subroutine read_wfc3ir_flt(FILENAME,pixr) implicit none real pixr(1014,1014) real pix(1014,1014), pmax integer*2 pox(1014,1014) integer pux(1014,1014) integer i, j integer ii,jj real xsum, ysum, psum, msky real xbar, ybar, mag real mbar_sky character*80 FILENAME character*70 INFO(10) common / fitsinfo / INFO real EXPTIME integer k real pval integer ival real pixi(1014,1014) real pixh(1014,1014) integer*2 pixo(1014,1014) integer*2 pixi2(1014,1014) real rsum real dd integer Ls real plist(9) integer NIT, NFIXN, NFIXT byte pixi2b(1014,1014,16), b integer*2 i2 do i = 0001, 1014 do j = 0001, 1014 pixi(i,j) = 0. pixi2(i,j) = 0 enddo enddo call readfits_r4e(FILENAME,pixi ,1014,1014,1) call readfits_i2e(FILENAME,pixi2,1014,1014,3) do i = 0001, 1014 do j = 0001, 1014 i2 = pixi2(i,j) pixi2b(i,j,16) = i2/32768 i2 = i2-i2/32768*32768 pixi2b(i,j,15) = i2/16384 i2 = i2-i2/16384*16384 pixi2b(i,j,14) = i2/08192 i2 = i2-i2/08192*08192 pixi2b(i,j,13) = i2/04096 i2 = i2-i2/04096*04096 pixi2b(i,j,12) = i2/02048 i2 = i2-i2/02048*02048 pixi2b(i,j,11) = i2/01024 i2 = i2-i2/01024*01024 pixi2b(i,j,10) = i2/00512 i2 = i2-i2/00512*00512 pixi2b(i,j,09) = i2/00256 i2 = i2-i2/00256*00256 pixi2b(i,j,08) = i2/00128 i2 = i2-i2/00128*00128 pixi2b(i,j,07) = i2/00064 i2 = i2-i2/00064*00064 pixi2b(i,j,06) = i2/00032 i2 = i2-i2/00032*00032 pixi2b(i,j,05) = i2/00016 i2 = i2-i2/00016*00016 pixi2b(i,j,04) = i2/00008 i2 = i2-i2/00008*00008 pixi2b(i,j,03) = i2/00004 i2 = i2-i2/00004*00004 pixi2b(i,j,02) = i2/00002 i2 = i2-i2/00002*00002 pixi2b(i,j,01) = i2 c if (pixi2(i,j).ne.000) pixi(i,j) = -500 if (pixi2(i,j).eq.256) pixi(i,j) = 1e10 pixr(i,j) = pixi(i,j) enddo enddo do i = 1, 1014 do j = 1, 1014 pixh(i,j) = pixr(i,j) enddo enddo do i = 1, 1014 do j = 1, 1014 rsum = 0. if (pixr(i,j).ge._HIFLAG_) then do ii = max(0001,i-5), min(1014,i+5) do jj = max(0001,j-5), min(1014,j+5) dd = sqrt(1.*(i-ii)**2+(j-jj)**2) rsum = rsum + 1.0/(1+dd)**2*pixh(ii,jj) enddo enddo pixr(i,j) = pixr(i,j) + rsum endif enddo enddo do i = 0001, 1014 do j = 0001, 1014 rsum = sqrt((i-359.)**2+(j-55.)**2) if (rsum.lt.25) pixr(i,j) = -751 enddo enddo return end c c c subroutine readfits_r4e(FILE,pix,NDIMX,NDIMY,NEXTENU) implicit none character*80 FILE integer NDIMX,NDIMY real pix(NDIMX,NDIMY) integer NEXTENU character*80 FILEU character*70 INFO(10) common / fitsinfo / INFO integer naxes integer laxis(3) common/laxis3_/laxis integer NXF integer NYF character*8 field character*40 stream integer nbyte0 integer nbyteE integer nbyte1 integer nbyte2 integer nbper integer i,ios, k integer j character*2880 buffc byte buffb(2880) equivalence (buffc,buffb) integer ifirst, i1, i2 integer np1, np2, npt integer nextend integer nread real bscale, bzero integer bitpix character*70 HDR(25) common/HDR/HDR logical DIAG data DIAG /.false./ integer NEND FILEU = FILE do i = 75,2,-1 if (FILE(i:i+4).eq.'.fits') FILEU = FILE(1:i+4) enddo if (DIAG) then print*,'enter readfits_r4...' print*,'FILE: ',FILE(1:60) endif open(10,file=FILE,status='old', . err=900,recl=2880,form='UNFORMATTED', . access='DIRECT') if (DIAG) print*,'...opened' naxes = -1 laxis(1) = 1 laxis(2) = 1 laxis(3) = 1 nextend = 0 do i = 1, 10 INFO(i) = ' ' enddo BSCALE = 1.0 BZERO = 0.0 NEND = 0 i = 0 nread = 0 100 continue i = i + 1 if (DIAG) print*,'READREC: ',i read(10,rec=i,iostat=ios) buffc do k = 0, 35, 1 if (DIAG) write(*,'(i4,1x,i4,1x,a80)') . i,k,buffc(k*80+1:k*80+80) field = buffc(k*80+01:k*80+08) stream = buffc(k*80+11:k*80+51) if (field.eq.'NAXIS ') read(stream,*) naxes if (field.eq.'NAXIS1 ') read(stream,*) laxis(1) if (field.eq.'NAXIS2 ') read(stream,*) laxis(2) if (field.eq.'NAXIS3 ') read(stream,*) laxis(3) if (field.eq.'NEXTEND ') read(stream,*) nextend if (field.eq.'BITPIX ') read(stream,*) bitpix if (field.eq.'BSCALE ') read(stream,*) bscale if (field.eq.'BZERO ') read(stream,*) bzero if (field.eq.'EXPTIME ') INFO(1) = stream if (field.eq.'FILTER ') INFO(2) = stream if (field.eq.'FILTNAM1') INFO(2) = stream if (field.eq.'FILENAME') INFO(3) = stream if (field.eq.'DATE-OBS') INFO(4) = stream if (field.eq.'TIME-OBS') INFO(5) = stream if (field.eq.'DEC_TARG') INFO(6) = stream if (field.eq.'RA_TARG ') INFO(7) = stream if (field.eq.'PA_V3 ') INFO(8) = stream if (field.eq.'PROPOSID') INFO(9) = stream if (field.eq.'CRPIX1 ') HDR(01) = stream if (field.eq.'CRPIX2 ') HDR(02) = stream if (field.eq.'CRVAL1 ') HDR(03) = stream if (field.eq.'CRVAL2 ') HDR(04) = stream if (field.eq.'CTYPE1 ') HDR(05) = stream if (field.eq.'CTYPE2 ') HDR(06) = stream if (field.eq.'CD1_1 ') HDR(07) = stream if (field.eq.'CD1_2 ') HDR(08) = stream if (field.eq.'CD2_1 ') HDR(09) = stream if (field.eq.'CD2_2 ') HDR(10) = stream if (field.eq.'ORIENTAT') HDR(11) = stream if (field.eq.'PA_APER ') HDR(12) = stream if (field.eq.'PA_V3 ') HDR(13) = stream if (field.eq.'DATE-OBS') HDR(14) = stream if (field.eq.'TIME-OBS') HDR(15) = stream if (field.eq.'EXPTIME ') HDR(16) = stream if (field.eq.'ROOTNAME') HDR(17) = stream if (field.eq.'TARGNAME') HDR(18) = stream if (field.eq.'RA_TARG ') HDR(19) = stream if (field.eq.'DEC_TARG') HDR(20) = stream if (field.eq.'PROPOSID') HDR(21) = stream if (field.eq.'FILTER ') HDR(22) = stream if (field.eq.'FILTER1 ') HDR(22) = stream if (field.eq.'FILTER2 ') HDR(23) = stream if (field.eq.'VAFACTOR') HDR(24) = stream if (field.eq.'END ') then !print*,' ---> NEND: ',NEND,NEXTENU if (NEND.eq.NEXTENU) goto 101 NEND = NEND + 1 endif enddo goto 100 101 continue nread = nread + 1 if (DIAG) then print*,' ' print*,'----------------------------------------' print*,' NREAD: ',nread print*,'NEXTEND: ',nextend print*,' NAXIS: ',naxes print*,' LAXIS: ',laxis(1),laxis(2),laxis(3) print*,' BITPIX: ',bitpix print*,' BSCALE: ',bscale print*,' BZERO: ',bzero print*,' NDIMX: ',NDIMX print*,' NDIMY: ',NDIMY print*,' ' endif ifirst = i+1 i1 = i i2 = i NXF = laxis(1) NYF = laxis(2) nbper = 4*laxis(1)*laxis(2) npt = laxis(1)*laxis(2) nbyte1 = 1 nbyte2 = nbper i1 = i+1 + nbyte1/2880 i2 = i+1 + nbyte2/2880 if (BITPIX.ne.-32) then print*,'prob' stop endif do i = i1, i2, 1 read(10,rec=i,iostat=ios) buffc nbyte0 = (i-ifirst)*2880+ 1 nbyteE = (i-ifirst)*2880+2880 np1 = (nbyte0-nbyte1)/4 + 1 np2 = (nbyteE-nbyte1)/4 + 1 c call buff2pix_r4(buffb,pix,np1,npt) call buff2pix_r4_edge(buffb,pix,np1,npt, . NDIMX,NDIMY,laxis(1),laxis(2)) if (DIAG) write(*,1115) i,np1,np2,npt 1115 format(1x,i8,1x,i10,1x,i10,1x,i10) enddo close(10) return 900 continue print*,'READFITS ERROR' stop end c------------------------------------------------------ c c subroutine buff2pix_r4_edge(buff,pix,n1,nt, . NXP,NYP,NXF,NYF) implicit none c character buff(2880) byte buff(2880) integer NXP,NYP real pix(NXP,NYP) integer n1,nt integer NXF,NYF real pbuff(720) common /sneaky/pbuff integer i integer npu, nbu integer NX, NY byte b(4) real r equivalence(r,b) do i = 1, 720 npu = n1+i-1 nbu = (i-1)*4 if (npu.ge.1.and.npu.le.nt) then NX = npu - (npu-1)/NXF*NXF NY = 1 + (npu-1)/NXF if (NX.le.NXP.and.NY.le.NYP) then if (.not.(_LINUX_)) then b(1) = buff(nbu+1) b(2) = buff(nbu+2) b(3) = buff(nbu+3) b(4) = buff(nbu+4) endif if ((_LINUX_)) then b(4) = buff(nbu+1) b(3) = buff(nbu+2) b(2) = buff(nbu+3) b(1) = buff(nbu+4) endif if (npu.ge.1.and.npu.le.nt) pix(NX,NY) = r endif endif enddo return end c---------------------------------------------------------- c c c subroutine readfits_i2e(FILE,pix,NDIMX,NDIMY,NEXTENU) implicit none character*80 FILE integer NDIMX,NDIMY integer NEXTENU integer*2 pix(NDIMX,NDIMY) character*80 FILEU character*70 INFO(10) common / fitsinfo / INFO integer naxes integer laxis(3) integer NXF integer NYF character*8 field character*40 stream integer nbyte0 integer nbyteE integer nbyte1 integer nbyte2 integer nbper integer i,ios, k integer j character*2880 buffc byte buffb(2880) equivalence (buffc,buffb) integer ifirst, i1, i2 integer np1, np2, npt integer nextend integer nread real bscale, bzero integer bitpix character*70 HDR(25) common/HDR/HDR logical DIAG data DIAG /.false./ integer NEND integer ii FILEU = FILE do i = 75,2,-1 if (FILE(i:i+4).eq.'.fits') then FILEU = FILE(1:i+4) do ii = i+5, 80 FILEU(ii:ii) = ' ' enddo endif enddo if (DIAG) then print*,'enter readfits_i2e...' print*,'FILEU: ',FILEU(1:60) endif open(10,file=FILEU,status='old', . err=900,recl=2880,form='UNFORMATTED', . access='DIRECT') if (DIAG) print*,'...opened' naxes = -1 laxis(1) = 1 laxis(2) = 1 laxis(3) = 1 nextend = 0 do i = 1, 10 INFO(i) = ' ' enddo BSCALE = 1.0 BZERO = 0.0 NEND = -1 i = 0 nread = 0 100 continue i = i + 1 if (DIAG) print*,'READREC: ',i read(10,rec=i,iostat=ios) buffc do k = 0, 35, 1 if (DIAG) write(*,'(i4,1x,i4,1x,a80)') . i,k,buffc(k*80+1:k*80+80) field = buffc(k*80+01:k*80+08) stream = buffc(k*80+11:k*80+51) if (field.eq.'NAXIS ') read(stream,*) naxes if (field.eq.'NAXIS1 ') read(stream,*) laxis(1) if (field.eq.'NAXIS2 ') read(stream,*) laxis(2) if (field.eq.'NAXIS3 ') read(stream,*) laxis(3) if (field.eq.'NEXTEND ') read(stream,*) nextend if (field.eq.'BITPIX ') read(stream,*) bitpix if (field.eq.'BSCALE ') read(stream,*) bscale if (field.eq.'BZERO ') read(stream,*) bzero if (field.eq.'EXPTIME ') INFO(1) = stream if (field.eq.'FILTNAM1') INFO(2) = stream if (field.eq.'FILENAME') INFO(3) = stream if (field.eq.'DATE-OBS') INFO(4) = stream if (field.eq.'TIME-OBS') INFO(5) = stream if (field.eq.'DEC_TARG') INFO(6) = stream if (field.eq.'RA_TARG ') INFO(7) = stream if (field.eq.'DEC_DEG ') INFO(6) = stream if (field.eq.'RA_DEG ') INFO(7) = stream if (field.eq.'PA_V3 ') INFO(8) = stream if (field.eq.'PROPOSID') INFO(9) = stream if (field.eq.'CRPIX1 ') HDR(01) = stream if (field.eq.'CRPIX2 ') HDR(02) = stream if (field.eq.'CRVAL1 ') HDR(03) = stream if (field.eq.'CRVAL2 ') HDR(04) = stream if (field.eq.'CTYPE1 ') HDR(05) = stream if (field.eq.'CTYPE2 ') HDR(06) = stream if (field.eq.'CD1_1 ') HDR(07) = stream if (field.eq.'CD1_2 ') HDR(08) = stream if (field.eq.'CD2_1 ') HDR(09) = stream if (field.eq.'CD2_2 ') HDR(10) = stream if (field.eq.'ORIENTAT') HDR(11) = stream if (field.eq.'PA_APER ') HDR(12) = stream if (field.eq.'PA_V3 ') HDR(13) = stream if (field.eq.'DATE-OBS') HDR(14) = stream if (field.eq.'TIME-OBS') HDR(15) = stream if (field.eq.'EXPTIME ') HDR(16) = stream if (field.eq.'ROOTNAME') HDR(17) = stream if (field.eq.'TARGNAME') HDR(18) = stream if (field.eq.'RA_TARG ') HDR(19) = stream if (field.eq.'DEC_TARG') HDR(20) = stream if (field.eq.'RA_DEG ') HDR(19) = stream if (field.eq.'DEC_DEG ') HDR(20) = stream if (field.eq.'PROPOSID') HDR(21) = stream if (field.eq.'FILTER1 ') HDR(22) = stream if (field.eq.'FILTER2 ') HDR(23) = stream if (field.eq.'VAFACTOR') HDR(24) = stream if (field.eq.'CCDGAIN ') HDR(25) = stream if (field.eq.'END ') then NEND = NEND + 1 if (NEND.ge.1) goto 101 endif enddo goto 100 101 continue nread = nread + 1 if (DIAG) then print*,' ' print*,'----------------------------------------' print*,' NREAD: ',nread print*,'NEXTEND: ',nextend print*,' NAXIS: ',naxes print*,' LAXIS: ',laxis(1),laxis(2),laxis(3) print*,' BITPIX: ',bitpix print*,' BSCALE: ',bscale print*,' BZERO: ',bzero print*,' NDIMX: ',NDIMX print*,' NDIMY: ',NDIMY print*,' ' endif ifirst = i+1 i1 = i i2 = i NXF = laxis(1) NYF = laxis(2) nbper = 2*laxis(1)*laxis(2) npt = laxis(1)*laxis(2) nbyte1 = 1 nbyte2 = nbper i1 = i+1 + nbyte1/2880 i2 = i+1 + nbyte2/2880 c print*,'NEND: ',NEND,NEXTENU,BITPIX if (NEND.ne.NEXTENU) goto 100 if (BITPIX.ne.16) then print*,'prob' stop endif do i = i1, i2, 1 read(10,rec=i,iostat=ios) buffc nbyte0 = (i-ifirst)*2880+ 1 nbyteE = (i-ifirst)*2880+2880 np1 = (nbyte0-nbyte1)/2 + 1 np2 = (nbyteE-nbyte1)/2 + 1 call buff2pix_i2e(buffb,pix(1,1),np1,npt) if (DIAG) write(*,1115) i,np1,np2,npt,i/laxis(1) 1115 format(1x,i8,1x,i10,1x,i10,1x,i10,1x,i6.6) enddo close(10) return 900 continue print*,'READFITS_I2E ERROR' stop end subroutine buff2pix_i2e(buff,pix,n1,nt) implicit none byte buff(2880) integer*2 pix(*) integer n1,nt byte b(2) integer ii equivalence(ii,b) integer i, npu, nbu do i = 1, 1440 npu = n1+i-1 nbu = (i-1)*2 if (.not.(_LINUX_)) then b(1) = buff(nbu+1) b(2) = buff(nbu+2) endif if ((_LINUX_)) then b(2) = buff(nbu+1) b(1) = buff(nbu+2) endif if (npu.ge.1.and.npu.le.nt) pix(npu) = ii enddo return end c---------------------------------------- c c this will perform the inverse c distortion correction... c subroutine wfc3ir_cg(xc,yc,xr,yr) implicit none real*8 xc,yc real*8 xr,yr real*8 xt, yt real*8 dx, dy, dd integer NIT real*8 xru, yru xr = xc yr = yc xru = xc yru = yc NIT = 0 1 continue NIT = NIT + 1 call wfc3ir_gcX(xru,yru,xt,yt) dx = xc-xt dy = yc-yt dd = sqrt(dx**2+dy**2) xru = xru + (xc-xt)*0.5 yru = yru + (yc-yt)*0.5 if (dd.gt.0.001.and.NIT.lt.99) goto 1 11 continue xr = xru yr = yru return end c----------------------------------------------------- c c this just writes a real*4 fits image c subroutine writfits_r4(FILE,pix,PXDIMX,PXDIMY) implicit none character*(*) FILE integer PXDIMX,PXDIMY real pix(PXDIMX,PXDIMY) integer nbyte0 integer nbyteE integer nbyte1 integer nbyte2 integer nbper integer i,ios character*2880 buffc byte buffb(2880) equivalence (buffb,buffc) integer ifirst, i1, i2 integer np1, np2, npt integer k character*80 FILEU character*70 HDR(25) common/HDR/HDR FILEU = FILE do i = 75,2,-1 if (FILE(i:i+4).eq.'.fits') FILEU = FILE(1:i+4) enddo open(10,file=FILEU,status='unknown', . err=900,recl=2880,form='UNFORMATTED', . access='DIRECT') write(buffc( 0*80+1: 1*80),'(''SIMPLE = T'')') write(buffc( 1*80+1: 2*80),'(''BITPIX = -32'')') write(buffc( 2*80+1: 3*80),'(''NAXIS = '',8x,i12)') 2 write(buffc( 3*80+1: 4*80),'(''NAXIS1 = '',8x,i12)') PXDIMX write(buffc( 4*80+1: 5*80),'(''NAXIS2 = '',8x,i12)') PXDIMY write(buffc( 5*80+1: 6*80),'(''DATATYPE= '',9a)') . " 'REAL*4' " write(buffc(07*80+1:08*80),'(''COMMENT '',a05)') ' ' write(buffc(08*80+1:09*80),'(''COMMENT '',a05)') ' ' write(buffc(09*80+1:10*80),'(''CRPIX1 = '',a70)') HDR(01) write(buffc(10*80+1:11*80),'(''CRPIX2 = '',a70)') HDR(02) write(buffc(11*80+1:12*80),'(''CRVAL1 = '',a70)') HDR(03) write(buffc(12*80+1:13*80),'(''CRVAL2 = '',a70)') HDR(04) write(buffc(13*80+1:14*80),'(''CTYPE1 = '',a70)') HDR(05) write(buffc(14*80+1:15*80),'(''CTYPE2 = '',a70)') HDR(06) write(buffc(15*80+1:16*80),'(''CD1_1 = '',a70)') HDR(07) write(buffc(16*80+1:17*80),'(''CD1_2 = '',a70)') HDR(08) write(buffc(17*80+1:18*80),'(''CD2_1 = '',a70)') HDR(09) write(buffc(18*80+1:19*80),'(''CD2_2 = '',a70)') HDR(10) write(buffc(19*80+1:20*80),'(''ORIENTAT= '',a70)') HDR(11) write(buffc(20*80+1:21*80),'(''PA_APER = '',a70)') HDR(12) write(buffc(21*80+1:22*80),'(''PA_V3 = '',a70)') HDR(13) write(buffc(22*80+1:23*80),'(''DATE-OBS= '',a70)') HDR(14) write(buffc(23*80+1:24*80),'(''TIME-OBS= '',a70)') HDR(15) write(buffc(24*80+1:25*80),'(''EXPTIME = '',a70)') HDR(16) write(buffc(25*80+1:26*80),'(''ROOTNAME= '',a70)') HDR(17) write(buffc(26*80+1:27*80),'(''TARGNAME= '',a70)') HDR(18) write(buffc(27*80+1:28*80),'(''RA_TARG = '',a70)') HDR(19) write(buffc(28*80+1:29*80),'(''DEC_TARG= '',a70)') HDR(20) write(buffc(29*80+1:30*80),'(''PROPOSID= '',a70)') HDR(21) write(buffc(30*80+1:31*80),'(''FILTER1 = '',a70)') HDR(22) write(buffc(31*80+1:32*80),'(''FILTER2 = '',a70)') HDR(23) write(buffc(33*80+1:34*80),'(''VAFACTOR= '',a70)') HDR(24) write(buffc(32*80+1:33*80),'(''CCDGAIN = '',a70)') HDR(25) write(buffc(33*80+1:34*80),'(''COMMENT '',a05)') ' ' write(buffc(34*80+1:35*80),'(''COMMENT '',a05)') ' ' write(buffc(35*80+1:36*80),'(''END '')') do k = 00, 34 if (buffc(k*80+11:k*80+30).eq.' ') . write(buffc(k*80+01:k*80+80),'(80('' ''))') enddo write(10,rec=i,iostat=ios) buffc ifirst = i+1 i1 = i i2 = i nbper = 4*PXDIMX*PXDIMY npt = PXDIMX*PXDIMY nbyte1 = 1 nbyte2 = nbper i1 = i+1 + nbyte1/2880 i2 = i+1 + nbyte2/2880 do i = i1, i2, 1 nbyte0 = (i-ifirst)*2880+ 1 nbyteE = (i-ifirst)*2880+2880 np1 = (nbyte0-nbyte1)/4 + 1 np2 = (nbyteE-nbyte1)/4 + 1 call pix2buff_r4(buffb,pix,np1,npt) write(10,rec=i,iostat=ios) buffc enddo close(10) return 900 continue print*,'writfits_r4.f ERROR' print*,' FILEU: ',FILEU stop end c------------------------------------------------------- c c subroutine pix2buff_r4(buff,pix,n1,nt) implicit none byte buff(2880) real*4 pix(*) integer n1,nt byte b(4) real*4 r equivalence(r,b) integer i, npu, nbu do i = 1, 720 npu = n1+i-1 nbu = (i-1)*4 if (npu.ge.1.and.npu.le.nt) r = pix(npu) if (.not.(_LINUX_)) then buff(nbu+1) = b(1) buff(nbu+2) = b(2) buff(nbu+3) = b(3) buff(nbu+4) = b(4) endif if ((_LINUX_)) then buff(nbu+1) = b(4) buff(nbu+2) = b(3) buff(nbu+3) = b(2) buff(nbu+4) = b(1) endif enddo return end