C+ C NAME: C program ashi_flat_field C PURPOSE: C To read in ASHI image grd files and provide a flat field from them. C C CALLING SEQUENCE: C ./ashi_flat_field C INPUTS: C C C FUNCTIONS/SUBROUTINES: C C MODIFICATION HISTORY: C September 2023, Bernard Jackson (UCSD) C- program ashi_flat_field parameter(NXY = 2048) ! Maximum number of x or y values in each grd image parameter(NOBS = 200) ! Potential number of valid images to flat field parameter(ihead = 5) ! # of lines in the GRD file header real aimagexy (NXY,NXY), ! Original image read in floating point if there is one & aimage (NXY,NXY,NOBS), ! images to combine & pcfix (NXY,NXY,NOBS), ! the fractional fix to each image & dfixi (NXY,NXY,NOBS), ! the image that is fixed - should be the same as the mean & dhmean (NXY,NXY,NOBS), ! the mean value of the +-5 box & aimagesm (NXY,NXY,NOBS), ! smoothed images that tests the images go smoother & pcfixsm (NXY,NXY,NOBS) ! the fractional fix to each smoothed image real pcfixxy (NXY,NXY), & pcfixxysm (NXY,NXY), & dhmeanxy (NXY,NXY), & dhmeanxysm(NXY,NXY), & tpcfixxy (NXY,NXY), & tpcfixxysm(NXY,NXY), & dtpcfixxy (NXY,NXY), & anhmean (NXY,NXY), & anhmeansm (NXY,NXY) integer imagexy (NXY,NXY) ! original image read in in ADU character cfile(NOBS)*52 ! Name of each input file to flat field character cfilel*66 /'/home/soft/mbracamontes/ASHI_Rot_data/Lab_Flats/Lab_Flat_00000.grd'/ ! Location and name of a file to flat field character cfilem*69 /'/home/soft/mbracamontes/ASHI_Rot_data/Lab_Flats/Mirror_Flat_00000.grd'/ character cfilei*69 character cfileoo*57 /'/home/bjackson/soft/for/ASHI_flt_files/ASHI_000_norm.grd'/ ! output file character cfileos*57 /'/home/bjackson/soft/for/ASHI_flt_files/ASHI_000_smoo.grd'/ ! output file character cfileom*57 /'/home/bjackson/soft/for/ASHI_flt_files/ASHI_000_mean.grd'/ ! output file character cfileon_m*57 /'/home/bjackson/soft/for/ASHI_flt_files/ASHI_000_n_mn.grd'/ ! output file character cfileos_m*57 /'/home/bjackson/soft/for/ASHI_flt_files/ASHI_000_s_mn.grd'/ ! output file character cfileo*57 character cimghead(ihead)*9 ! Grd file 5 line header distcx = 1126.5 ! center of the image in X as determined by Matthew distcy = 1000.5 ! center of the image in Y as determined by Matthew dmax = 920.0 ! no sense in going too far beyond the image center dmin = 500.0 ! no sense in going too to close to the image cennter for the itype = 2 ! read in a grd file with <= 5 integer digits; output a file with F8.1 floating values NOBB = NOBS C Lab = 0 ! if Lab = 0, this is a Mirror flat Lab = 1 ! if Lab = 1, this is a Lab flat C Lao = 1 Lao = 2 ! output a F9.2 file NOBB = 111 NOBE = 191 do K=NOBB,NOBE do J=1,NXY do I=1,NXY aimage (I,J,K) = 0.0 pcfix (I,J,K) = 0.0 dfixi (I,J,K) = 0.0 dhmean (I,J,K) = 0.0 aimagesm (I,J,K) = 0.0 pcfixsm (I,J,K) = 0.0 if(K.eq.NOBB) then aimagexy (I,J) = 0.0 pcfixxy (I,J) = 0.0 pcfixxysm (I,J) = 0.0 dhmeanxy (I,J) = 0.0 dhmeanxysm(I,J) = 0.0 tpcfixxy (I,J) = 0.0 tpcfixxysm(I,J) = 0.0 dtpcfixxy (I,J) = 0.0 anhmean (I,J) = 0.0 anhmeansm (I,J) = 0.0 imagexy (I,J) = 0 end if end do end do end do do K=NOBB,NOBE inum1 = K if(Lab.eq.1.and.Lao.gt.0) then cfilei = cfilel if(inum1.ge.0 .and.inum1.le.9) write(cfilei(62:62),'(I1)') inum1 if(inum1.ge.10 .and.inum1.le.99) write(cfilei(61:62),'(I2)') inum1 if(inum1.ge.100.and.inum1.le.999) write(cfilei(60:62),'(I3)') inum1 end if if(Lab.eq.0.and.Lao.gt.0) then cfilei = cfilem if(inum1.ge.0 .and.inum1.le.9) write(cfilei(65:65),'(I1)') inum1 if(inum1.ge.10 .and.inum1.le.99) write(cfilei(64:65),'(I2)') inum1 if(inum1.ge.100.and.inum1.le.999) write(cfilei(63:65),'(I3)') inum1 end if write(*,'(A,1X,A)') ' cfilei',cfilei call read_grd_image(itype,cfilei,cfileo,NXY,dmax,cimghead,aimagexy,imagexy) do J=1,NXY do I=1,NXY aimage(I,J,K) = aimagexy(I,J) end do end do end do NV = 5 ! The box is plus and minus this limit in size alol = 5000.0 ! The lower ADU limit of the pixels accepted - about 6000.0 aupl = 15000.0 ! The upper ADU limit of the pixels accepted - about 14000.0 ped = 115.0 ! pedestal + dark current (assumed constant) tfixxy = 0.0 antfixxy = 0.0 do K=NOBB,NOBE do J=1,NXY do I=1,NXY xdis = distcx-float(I) ydis = distcy-float(J) rdis = sqrt(xdis**2 + ydis**2) if(rdis.le.dmax.and.rdis.ge.dmin) then dhum = 0.0 anum = 0.0 do JJ=J-NV,J+NV do II=I-NV,I+NV if(JJ.ge.1.and.JJ.le.NXY.and.II.ge.1.and.II.le.NXY) then if(aimage(II,JJ,K).gt.alol.and.aimage(II,JJ,K).lt.aupl) then dhum = dhum + aimage(II,JJ,K) anum = anum + 1.0 end if end if end do end do if(anum.ne.0) then dhmean(I,J,K) = (dhum/anum)-ped ! mean value of the +-5 box pcfix (I,J,K) = (aimage(I,J,K)-dhmean(I,J,K))/dhmean(I,J,K) ! fractional difference of the current pixel I, J, K from the mean C How to change this one image so it is fixed using the fractional mean dfixi (I,J,K) = aimage(I,J,K) - pcfix(I,J,K)*dhmean(I,J,K) ! the fixed image that is now smooth C pcfixxy (I,J) = pcfixxy (I,J) + pcfix(I,J,K)*dhmean(I,J,K) ! sum all the image fixes times the mean ! Oops C pcfixxy (I,J) = pcfixxy (I,J) - pcfix(I,J,K)*dhmean(I,J,K) ! sum all the image fixes times the mean pcfixxy (I,J) = pcfixxy (I,J) - pcfix(I,J,K) ! sum all the image pixel fixes anhmean(I,J) = anhmean(I,J) + 1.0 ! find the number of all the summed image means III = 500 JJJ = 500 if(I.eq.III.and.J.eq.JJJ) & write(*,'(I3,F9.3,F10.7,F9.3,3F10.2)') & K,dhmean(III,JJJ,K),pcfix(III,JJJ,K),dfixi(III,JJJ,K),pcfixxy(III,JJJ),anhmean(III,JJJ),anum end if end if end do end do end do C Now zero sections of the image that are fixed, but have blemishes C Set 1 1-99 image grid locations to zero if(NOBB.lt.100) then thea1x = 704.0 thea1y = 1511.0 sizet1 = 1.0 thea2x = 494.0 thea2y = 851.0 sizet2 = 5.0 thea3x = 1574.0 thea3y = 549.0 sizet3 = 1.0 thea4x = 1746.0 thea4y = 1203.0 sizet4 = 5.0 rspt1x = 724.0 rspt1y = 1708.7 rdiam1 = 15.0 rspt2x = 1723.4 rspt2y = 1429.2 rdiam2 = 19.0 xcheck = 1160.0 ycheck = 1730.0 sizet5 = 75.0 end if C Set 2 1-99 image grid locations to zero if(NOBB.ge.101) then thea1x = 1044.5 thea1y = 1640.0 sizet1 = 1.0 thea2x = 292.5 thea2y = 1294.0 sizet2 = 5.0 thea3x = 1240.4 thea3y = 377.5 sizet3 = 1.0 thea4x = 1964.4 thea4y = 788.5 sizet4 = 5.0 rspt1x = 724.0 rspt1y = 1708.7 rdiam1 = 15.0 rspt2x = 1723.4 rspt2y = 1429.2 rdiam2 = 19.0 xcheck = 1160.0 ycheck = 1730.0 sizet5 = 75.0 end if call pixelcon(2,rspt1x,rspt1y,thea4x,thea4y,size,rdiam1,dhmean,NXY,NOBB,NOBE,distcx,distcy) ! eliminate a large bad spot call pixelcon(2,rspt2x,rspt2y,thea4x,thea4y,size,rdiam2,dhmean,NXY,NOBB,NOBE,distcx,distcy) ! eliminate a large bad spot call pixelcon(3,x,y,thea1x,thea1y,sizet1,rdiam,dhmean,NXY,NOBB,NOBE,distcx,distcy) ! eliminate a bad sector call pixelcon(3,x,y,thea2x,thea2y,sizet2,rdiam,dhmean,NXY,NOBB,NOBE,distcx,distcy) ! eliminate a bad sector call pixelcon(3,x,y,thea3x,thea3y,sizet3,rdiam,dhmean,NXY,NOBB,NOBE,distcx,distcy) ! eliminate a bad sector call pixelcon(3,x,y,thea4x,thea4y,sizet4,rdiam,dhmean,NXY,NOBB,NOBE,distcx,distcy) ! eliminate a bad sector call pixelcon(4,xcheck,ycheck,thea4x,thea4y,sizet5,rdiam,dhmean,NXY,NOBB,NOBE,distcx,distcy) ! check a smooth location call pixelcon(4,xcheck,ycheck,thea4x,thea4y,sizet5,rdiam,aimage,NXY,NOBB,NOBE,distcx,distcy) ! provide an unsmoothed image location C We think each image pixel fractional fix is going to be the same +- value relative to the mean at that portion of the image C Therefore sum of all the pixel fractional fixes at that image location should be divided by the total mean values you found for the correction do J=1,NXY do I=1,NXY if(anhmean(I,J).ne.0.0) then tpcfixxy (I,J) = pcfixxy (I,J)/anhmean(I,J) ! these are the mean fractional fixes of all of the individual NOB pixels end if end do end do III = 500 JJJ = 500 print*, 'sample of all fixes', pcfixxy(III,JJJ), anhmean(III,JJJ), tpcfixxy(III,JJJ) C Now that we have the mean fractional fix tpcfixxy for each pixel, we can see how much each of these fixes varies from each other tfix = 0.0 antfix = 0.0 do K=NOBB,NOBE do J=1,NXY do I=1,NXY if(dhmeanxy(I,J).ne.0.0) then tfix = tfix + (pcfix(i,k,K) - tpcfixxy(I,J)) antfix = antfix + 1.0 end if end do end do end do print*, 'tfix,antfix',tfix,antfix tfix = tfix/antfix ! the average fractional fix difference from the mean fix - ideally zero write(*,'(A,I4,A,F8.5)') ' The average fractional fix of ',NOBE-NOBB,' images is ',tfix C Now fix all the NOB images by the mean fix of each pixel do K=NOBB,NOBE do J=1,NXY do I=1,NXY if(dhmean(I,J,K).ne.0.0) then aimagesm(I,J,K) = aimage(I,J,K) - tpcfixxy(I,J)*dhmean(I,J,K) ! the fixed image that is now more smooth end if end do end do end do C Now see if the variations of each smoothed pixel on the average is smaller do K=NOBB,NOBE do J=1,NXY do I=1,NXY if(dhmean(I,J,K).ne.0.0) then pcfixsm (I,J,K) = (aimagesm(I,J,K)-dhmean(I,J,K))/dhmean(I,J,K) ! fractional difference of the smoothed pixel I, J, K from the mean C pcfixxysm (I,J) = pcfixxysm (I,J) + pcfixsm(I,J,K)*dhmean(I,J,K) ! sum all the image fixes times the mean C dhmeanxysm(I,J) = dhmeanxysm(I,J) + dhmean(I,J,K) ! sum all the image means pcfixxysm (I,J) = pcfixxysm (I,J) + pcfixsm(I,J,K) ! sum all the image fixes anhmeansm(I,J) = anhmeansm (I,J) + 1.0 ! find the number of all the sumed image means III = 500 JJJ = 500 if(I.eq.III.and.J.eq.JJJ) & write(*,'(I3,F9.3,F10.7,9X,2F10.2)'), & K,dhmean(III,JJJ,K),pcfixsm(III,JJJ,K), pcfixxysm(III,JJJ),anhmeansm(III,JJJ) end if end do end do end do do J=1,NXY do I=1,NXY if(dhmeanxysm(I,J).ne.0.0) then tpcfixxysm (I,J) = pcfixxysm(I,J)/dhmeanxysm(I,J) ! these are the mean fractional fixes of all of the smoothed individual NOB pixels tpcfixxysm (I,J) = pcfixxysm(I,J)/anhmeansm(I,J) ! these are the mean fractional fixes of all of the smoothed individual NOB pixels end if end do end do III = 500 JJJ = 500 C print*, 'sample of all smoothed fixes', pcfixxysm(III,JJJ), dhmeanxysm(III,JJJ), tpcfixxysm(III,JJJ) print*, 'sample of all smoothed fixes', pcfixxysm(III,JJJ), anhmeansm(III,JJJ), tpcfixxysm(III,JJJ) C Now that we have the mean fractional fix tpcfixxysm for each smoothed pixel, we can see how much each of these fixes varies from each other tfixsm = 0.0 antfixsm = 0.0 do K=NOBB,NOBE do J=1,NXY do I=1,NXY C if(dhmeanxysm(I,J).ne.0.0) then if(anhmeansm(I,J).ne.0.0) then tfixsm = tfixsm + (pcfixsm(I,J,K) - tpcfixxysm(I,J)) antfixsm = antfixsm + 1.0 end if end do end do end do print*, 'tfixsm,antfixsm',tfixsm,antfixsm tfixsm = tfixsm/antfixsm ! the average fractional fix difference from the mean fix - ideally zero write(*,'(A,I4,A,F10.8)') ' The average fractional fix of ',NOBE-NOBB,' smoothed images is ',tfixsm print*,' ' C OK, so now write out an image that is not smoothed, and the same one that is smoothed, and check to see if it is better C Do this with an image in the center of the stack used to make the smoothed image NOBI = (NOBE-NOBB)/2 + NOBB - 1 print*, 'The image to check is ',NOBI C Sample Original Image do J=1,NXY do I=1,NXY aimagexy(I,J) = 0.0 aimagexy(I,J) = aimage(I,J,NOBI) + ped end do end do inum1 = NOBI dmax = 3000. if(Lao.eq.2) then cfileo = cfileoo if(inum1.ge.0 .and.inum1.le.9) write(cfileo(47:47),'(I1)') inum1 if(inum1.ge.10 .and.inum1.le.99) write(cfileo(46:47),'(I2)') inum1 if(inum1.ge.100.and.inum1.le.999) write(cfileo(45:47),'(I3)') inum1 end if write(*,'(A,1X,A)') ' cfileoo',cfileo call read_grd_image(15,cfilei,cfileo,NXY,dmax,cimghead,aimagexy,imagexy) C Sample Smoothed image do J=1,NXY do I=1,NXY aimagexy(I,J) = 0.0 aimagexy(I,J) = aimagesm(I,J,NOBI) + ped end do end do if(Lao.eq.2) then cfileo = cfileos if(inum1.ge.0 .and.inum1.le.9) write(cfileo(47:47),'(I1)') inum1 if(inum1.ge.10 .and.inum1.le.99) write(cfileo(46:47),'(I2)') inum1 if(inum1.ge.100.and.inum1.le.999) write(cfileo(45:47),'(I3)') inum1 end if write(*,'(A,1X,A)') ' cfileos',cfileo call read_grd_image(15,cfilei,cfileo,NXY,dmax,cimghead,aimagexy,imagexy) C Sample Mean image do J=1,NXY do I=1,NXY aimagexy(I,J) = 0.0 aimagexy(I,J) = dhmean(I,J,NOBI) + ped end do end do if(Lao.eq.2) then cfileo = cfileom if(inum1.ge.0 .and.inum1.le.9) write(cfileo(47:47),'(I1)') inum1 if(inum1.ge.10 .and.inum1.le.99) write(cfileo(46:47),'(I2)') inum1 if(inum1.ge.100.and.inum1.le.999) write(cfileo(45:47),'(I3)') inum1 end if write(*,'(A,1X,A)') ' cfileom',cfileo call read_grd_image(15,cfilei,cfileo,NXY,dmax,cimghead,aimagexy,imagexy) C Sample Normal minus Mean image do J=1,NXY do I=1,NXY aimagexy(I,J) = 0.0 aimagexy(I,J) = aimage(I,J,NOBI) - dhmean(I,J,NOBI) end do end do if(Lao.eq.2) then cfileo = cfileon_m if(inum1.ge.0 .and.inum1.le.9) write(cfileo(47:47),'(I1)') inum1 if(inum1.ge.10 .and.inum1.le.99) write(cfileo(46:47),'(I2)') inum1 if(inum1.ge.100.and.inum1.le.999) write(cfileo(45:47),'(I3)') inum1 end if write(*,'(A,1X,A)') ' cfileon_m',cfileo call read_grd_image(15,cfilei,cfileo,NXY,dmax,cimghead,aimagexy,imagexy) C Sample Smooth minus Mean image do J=1,NXY do I=1,NXY aimagexy(I,J) = 0.0 aimagexy(I,J) = aimagesm(I,J,NOBI) - dhmean(I,J,NOBI) end do end do if(Lao.eq.2) then cfileo = cfileos_m if(inum1.ge.0 .and.inum1.le.9) write(cfileo(47:47),'(I1)') inum1 if(inum1.ge.10 .and.inum1.le.99) write(cfileo(46:47),'(I2)') inum1 if(inum1.ge.100.and.inum1.le.999) write(cfileo(45:47),'(I3)') inum1 end if write(*,'(A,1X,A)') ' cfileos_m',cfileo call read_grd_image(15,cfilei,cfileo,NXY,dmax,cimghead,aimagexy,imagexy) stop end