PROGRAM TIFTOFITS c rku USE AVDEF c use tiffwrite c use tiffread c end rku C C CONVERTS A TIFF FILE TO FITS. BLACK & WHITE 16 BIT TRANSPARENCY C SCAN FROM MICROTEK/PHOTOSHOP. C COMMON/EPHE/P,B0,DEC0,SINCLN,COSCLN,SINB0,COSB0,SINP,COSP,PHI,FE, 1 H0,CLONG0, ZD CHARACTER*80 MTIF, MLIN, MFIT, MDUSTFIT, MSLOPEFIT, MRMSFIT, mtxt 1 ,mprodfit character*80 mfile_list,mheader, mlin0,m_inpr,m_outpr character*30 myearcase,mmonthcase,mrootin,mrootout,mrootdiag 1 ,mrootprod character*30 myrmocase character*60 mrtin,mrtout,mrtdiag,mrtprod integer*2 narg1,narg2 CHARACTER*4 MRESU,myr CHARACTER*1 MYON character*9 daterun character*10 dateid character*10 timerun character*5 timeid character*40 caseid character*70 mtifkwd_v character*20 mtifiden,mtifiden_in,mtif_integer,mtifkwd_v0 character*20 mtifid(38) character*25 mtifid0 character mtifkwd(38)*8,mtifkwd_value(38)*70 character*2 xps,xms,yps,yms logical allocated,kwd_is_string(38) logical carry_on,grad_test,rmin_test,rlim_test,dust_test, 1 blank_test,line_test,edge_test logical grad_test0,rmin_test0,rlim_test0,dust_test0, 1 blank_test0,line_test0,edge_test0 INTEGER*2 IBUF(80000), IREF(2), I2MAD(800), IBFO(1440), NUL, 1 IPIC(100000000) DIMENSION ITAGTAB(12), IMAD(400), BUFO(720), AR(5), ARI(40000), 1 ARO(40000), WEI(40000),fpic(10000000) EQUIVALENCE (IDATA,IREF), (IMAD,I2MAD), (BUFO,IBFO) INTEGER*4 I4DAR(5000) INTEGER*2 I2DAR(10000) BYTE I1DAR(20000) DIMENSION RDAR(5000) EQUIVALENCE (I4DAR,I2DAR), (I4DAR,I1DAR), (I4DAR,RDAR) DATA ITAGTAB/ 256, 257, 258, 259, 262, 270, 273, 277, 278, 1 282, 283, 296/ DATA NUL/ Z8000/ data rzero/1000.0/ data rmin/0.95/ data ams_fac/0.6/ data edge_fac/1.0/ real line_veto data line_veto/0.3/ data dust_veto/10.0/ data scattered_light_factor/0.7/ data grad_factor/0.2/ data offset_distance_minimum/0.15/ data offset_distance_maximum/0.45/ data scmax/1.5/ data num_input_tif_items/38/ data mtifid/ 1'Year', 1'Month', 1'Day', 1'Image_No', 1'ScanId', 1'plate_no', 1'PlateTypeCode', 1'ObsType', 1'time', 1'Image_Qual', 1'Logbook_Page', 1'LogComments', 1'TagComments', 1'ExtRes', 1'ExtQual', 1'FileName', 1'ImageWidth', 1'ImageLength', 1'ExtSubSample', 1'ScanResUnit', 1'ScanXres', 1'ScanYres', 1'PosInScan', 1'HorPosNum', 1'HorPos', 1'PlateXoffset', 1'PlateYoffset', 1'ImageXoffset', 1'ImageYoffset', 1'LogDataEntryDate', 1'LogSerial', 1'TagDataEntryDate', 1'TagSerial', 1'ExtFileLoc', 1'ScanFileLoc', 1'ScanFrame', 1'ExtDate', 1'Serial'/ data kwd_is_string/ 1 .FALSE., !'Year', 1 .FALSE., !'Month', 1 .FALSE., !'Day', 1 .FALSE., !'Image_No', 1 .TRUE., !'ScanId', 1 .TRUE., !'plate_no', 1 .FALSE., !'PlateTypeCode', 1 .TRUE., !'ObsType', 1 .TRUE., !'time', 1 .FALSE., !'Image_Qual', 1 .TRUE., !'Logbook_Page', 1 .TRUE., !'LogComments', 1 .TRUE., !'TagComments', 1 .FALSE., !'ExtRes', 1 .FALSE., !'ExtQual', 1 .TRUE., !'FileName', 1 .FALSE., !'ImageWidth', 1 .FALSE., !'ImageLength', 1 .FALSE., !'ExtSubSample', 1 .FALSE., !'ScanResUnit', 1 .FALSE., !'ScanXresRaw', 1 .FALSE., !'ScanYresRaw', 1 .FALSE., !'PosInScan', 1 .FALSE., !'HorPosNum', 1 .TRUE., !'HorPos C', 1 .FALSE., !'PlateXoffset', 1 .FALSE., !'PlateYoffset', 1 .FALSE., !'ImageXoffset', 1 .FALSE., !'ImageYoffset', 1 .TRUE., !'LogDataEntryDate', 1 .FALSE., !'LogSerial', 1 .TRUE., !'TagDataEntryDate', 1 .FALSE., !'TagSerial', 1 .TRUE., !'ExtFileLoc', 1 .TRUE., !'ScanFileLoc', 1 .FALSE., !'ScanFrame', 1 .TRUE., !'ExtDate', 1 .FALSE./ !'Serial'/ data mtifkwd/ 1'YEAR ', 1'MONTH ', 1'DAY ', 1'IMAGE_NO', 1'SCANID ', 1'PLATE_NO', 1'PLTYPECD', 1'OBSTYPE ', 1'TIME ', 1'IMAGQUAL', 1'LGBOOKPG', 1'LGCOMMNT', 1'TGCOMMNT', 1'EXTRES ', 1'EXTQUAL ', 1'FILENAME', 1'IMAGWDTH', 1'IMAGLNGT', 1'EXSUBSAM', 1'SCRESUNI', 1'SCANXRES', 1'SCANYRES', 1'POSINSCN', 1'HORPOSNM', 1'HORPOS ', 1'PLXOFFST', 1'PLYOFFST', 1'IMXOFFST', 1'IMYOFFST', 1'LGDTDATE', 1'LGSERIAL', 1'TGDTDATE', 1'TGSERIAL', 1'EXFILLOC', 1'SCFILLOC', 1'SCFRAME ', 1'EXTDATE ', 1'SERIAL '/ c rku INTEGER*2, ALLOCATABLE::ARR(:,:) , avearr(:,:) , difarr(:,:) , 1 mask(:,:), arr0(:,:) , slopeave(:) parameter nptave=20, nptrad=32, nptms=25, nptrow=35, nptdiag = 30 logical good, ijswap, ijdiagonal, reject type pr integer ii,jj end type pr type (pr), dimension(nptave):: ij data (ij(k),k=1,nptave) x /pr(-2,-2), pr(-2,-1),pr(-2,0) , pr(-2,1), pr(-2,2), x pr(-1,2),pr(0,2),pr(1,2),pr(2,2),pr(2,1), x pr(2,0),pr(2,-1),pr(2,-2),pr(1,-2),pr(0,-2), 2 pr(-1,-2),pr(-1,1),pr(1,1),pr(1,-1),pr(-1,-1)/ type (pr), dimension(nptrad):: ijr data (ijr(k) , k = 1, nptrad) x /pr(-3,-3),pr(-4,-3),pr(-4,-2),pr(-4,-1),pr(-4,0), x pr(-4,1),pr(-4,2),pr(-4,3),pr(-3,3),pr(-3,4),pr(-2,4), x pr(-1,4),pr(0,4),pr(1,4),pr(2,4),pr(3,4),pr(3,3), x pr(4,3),pr(4,2),pr(4,1),pr(4,0),pr(4,-1),pr(4,-2), x pr(4,-3),pr(3,-3),pr(3,-4),pr(2,-4),pr(1,-4),pr(0,-4), x pr(-1,-4),pr(-2,-4),pr(-3,-4)/ data testfac /4.0/ type (pr), dimension(nptms):: ijms data (ijms(k) , k = 1 , nptms) x /pr(-2,-2),pr(-2,-1),pr(-2,0),pr(-2,1),pr(-2,2), x pr(-1,-2),pr(-1,-1),pr(-1,0),pr(-1,1),pr(-1,2), x pr(0,-2),pr(0,-1), pr(0,0), pr(0,1),pr(0,2), x pr(1,-2),pr(1,-1), pr(1,0), pr(1,1),pr(1,2), x pr(2,-2),pr(2,-1), pr(2,0), pr(2,1),pr(2,2)/ type (pr), dimension(nptrow):: ijrow data (ijrow(k) , k = 1 , nptrow) x / x pr(1,-2),pr(2,-2),pr(3,-2),pr(4,-2),pr(5,-2),pr(6,-2),pr(7,-2), x pr(1,-1),pr(2,-1),pr(3,-1),pr(4,-1),pr(5,-1),pr(6,-1),pr(7,-1), x pr(1,0),pr(2,0),pr(3,0),pr(4,0),pr(5,0),pr(6,0),pr(7,0), x pr(1, 1),pr(2, 1),pr(3, 1),pr(4, 1),pr(5, 1),pr(6,1),pr(7,1), x pr(1, 2),pr(2, 2),pr(3, 2),pr(4, 2),pr(5, 2),pr(6,2),pr(7,2) x / type (pr), dimension(nptdiag):: ijdiag data (ijdiag(k) , k = 1 , nptdiag) x / x pr(4,1),pr(5,2),pr(6,3), x pr(3,1),pr(4,2),pr(5,3),pr(6,4), x pr(2,1),pr(3,2),pr(4,3),pr(5,4),pr(6,5), x pr(1,1),pr(2,2),pr(3,3),pr(4,4),pr(5,5),pr(6,6), x pr(1,2),pr(2,3),pr(3,4),pr(4,5),pr(5,6), x pr(1,3),pr(2,4),pr(3,5),pr(4,6), x pr(1,4),pr(2,5),pr(3,6) x / allocatable xave(:), yave(:),xslopeave(:),yslopeave(:) c end rku c TYPE *, 'FILENAME for File List:' c ACCEPT 30, mfile_list 30 FORMAT(A) narg1=1 narg2=2 call getarg(narg1,myearcase) call getarg(narg2,mmonthcase) call getenv('inroot',mrootin) call getenv('outroot',mrootout) call getenv('diag',mrootdiag) call getenv('prodroot',mrootprod) ipos1=len_trim(mrootdiag)+1 ilen2=len_trim(myearcase) myrmocase=myearcase myrmocase(ilen2+1:ilen2+1)= '\' myrmocase(ilen2+2:) = mmonthcase ilen4= len_trim(myrmocase) ipos2=ipos1+ilen4 mfile_list(ipos2:)='\files.txt' mrtin=mrootin(1:len_trim(mrootin))//myrmocase(1:ilen4)//'\' mrtout=mrootout(1:len_trim(mrootout))//myrmocase(1:ilen4)//'\' mrtdiag=mrootdiag(1:len_trim(mrootdiag))//myrmocase(1:ilen4)//'\' mrtprod=mrootprod(1:len_trim(mrootprod))//myrmocase(1:ilen4)//'\' ilin=len_trim(mrtin) ilout=len_trim(mrtout) ildiag=len_trim(mrtdiag) ilprod=len_trim(mrtprod) mfile_list=mrtdiag(1:ildiag)//'files.txt' open(unit = 10, file = mfile_list,carriagecontrol = 'list', 1 status = 'old') allocated = .false. open(unit=11,file=mrtout(1:ilout)//'summary.txt',access='append', 1 carriagecontrol='list') c open(unit=12,file='header.txt',status='replace', c 1 carriagecontrol='list') 35 continue read(10,30,end=880) mtif type * , mtif if (len_trim(mtif).lt.5) go to 880 c 32 TYPE *, 'OBSERVATION DATE & TIME (MM/DD/YYYY, HH:MM):' c CALL RENOSB(5, AR, NNR, 5, '/, :') c IF(NNR.NE.5) GO TO 32 c IMO = AR(1) c IDA = AR(2) c IYR = AR(3) c TAV = AR(4) + AR(5)/60. c IF(IYR.LT.1800) GO TO 32 c IF(TAV.LT.14.) TAV = TAV + 8.0 c TYPE *, 'WAVELENGTH (=3933.682 CA, 0=DIRECT):' c CALL RENOSB(5, WAVEL, NNR, 1, ' ,') c IF(NNR.EQ.0) WAVEL = 3933.682 c TYPE *, 'IS ORIGINAL A NEGATIVE IMAGE (=Y)?' c ACCEPT 34, NC, MYON c 34 FORMAT(Q,A) IFLIP = 1 c IF(NC.EQ.0.OR.MYON.EQ.'Y'.OR.MYON.EQ.'y') IFLIP = 1 c TYPE *, 'INCREASE X AXIS LENGTH BY 6/5 (INTERPOLATE) TO', c 1 ' ROUNDIFY (=N)?' c ACCEPT 34, NC, MYON IASPEC = 0 c IF(NC.NE.0.AND.(MYON.EQ.'Y'.OR.MYON.EQ.'y')) THEN c IASPEC = 1 c IREBF = 1 c ENDIF OPEN(UNIT=1, FILE=mrtin(1:ilin)//MTIF, STATUS='OLD', 1 FORM='UNFORMATTED', 1 ACCESS='DIRECT', RECL=128) NC = INDEX(MTIF, '.') MFIT(1:) = mrtout(1:ilout)//MTIF(1:NC)//'fits' OPEN(UNIT=3, FILE=MFIT, STATUS='REPLACE', FORM='UNFORMATTED', 1 ACCESS='DIRECT', RECL=720) mdustfit(1:) = mrtdiag(1:ildiag)//mtif(1:nc-1)//'dust.fits' open(unit=4, file=mdustfit,status='REPLACE', form = 'unformatted', 1 access='direct', recl=720) mrmsfit(1:) = mrtdiag(1:ildiag)//mtif(1:nc-1)//'rms.fits' open(unit=2, file=mrmsfit,status='REPLACE', form = 'unformatted', 1 access='direct', recl=720) mslopefit(1:) = mrtdiag(1:ildiag)//mtif(1:nc-1)//'slope.fits' open(unit=7, file=mslopefit,status='REPLACE', 1 form = 'unformatted', 1 access='direct', recl=720) mprodfit(1:) = mrtprod(1:ilprod)//mtif(1:nc-1)//'.fits' open(unit=13, file=mprodfit,status='REPLACE', 1 form = 'unformatted', 1 access='direct', recl=720) mtxt(1:) = mrtdiag(1:ildiag)//mtif(1:nc)//'txt' OPEN(UNIT=8, FILE=mtxt, STATUS='REPLACE', 1 CARRIAGECONTROL='LIST') write(8,*) mtxt c mtxt(1:) = mtif(1:nc-1)//'slopes.txt' c OPEN(UNIT=9, FILE=mtxt, STATUS='REPLACE', c 1 CARRIAGECONTROL='LIST') c write(9,*) mtxt call getdat(iyearrun,imonthrun,idayrun) call gettim(ihourrun,iminrun,isecrun,ihundrun) call date_and_time(daterun,timerun) caseid = mtif(1:nc-1) dateid(1:4) =daterun(1:4) dateid(5:) = '-' dateid(6:7) = daterun(5:6) dateid(8:) = '-' dateid(9:10) = daterun(7:8) timeid(1:2) =timerun(1:2) timeid(3:) = ':' timeid(4:5) = timerun(3:4) caseid(25:34) = dateid caseid(36:40) = timeid c write(11,*) mtif(1:nc-1),' ',daterun,timerun CALL GETIFEL(0, IBUF, 10, 0, 1) c TYPE *, (IBUF(I),I=1,5) IF(IBUF(1).NE.18761) GO TO 40 ISWBF = 0 GO TO 50 40 IF(IBUF(1).NE.19789) GO TO 950 ISWBF = 1 CALL GETIFEL(0, IBUF, 10, ISWBF, 1) 50 IF(IBUF(2).NE.42) GO TO 950 IFDA = I2TOI4(IBUF(3)) CALL GETIFEL(IFDA, BUF, 2, ISWBF, 1) NIFD = IBUF(1) IFDA = IFDA + 2 c TYPE *, 'FILE DIR ADDR, NO. ENTRIES', IFDA, NIFD DO 320 I=1,NIFD CALL IFDNTRY(IFDA, ISWBF, 1, ITAG, ITYP, ICOUNT, IDATA, RDATA, 1 IVAF) c TYPE *, ITAG, ITYP, ICOUNT, IDATA, RDATA, IVAF DO 60 J=1,12 IF(ITAG.NE.ITAGTAB(J)) GO TO 60 NDX = J GO TO 70 60 CONTINUE GO TO 320 C TAG > 256 257 258 259 262 270 273 277 278 282 283 296 70 GO TO (80, 90, 100, 110, 120, 125, 130, 200, 210, 220, 230, 240) 1 , NDX C NO. PIX IN X. 80 NAXIS1 = IDATA GO TO 320 C NO. PIX IN Y. 90 NAXIS2 = IDATA GO TO 320 C NO. BITS PER PIXEL. 100 IBITPIX = IDATA IF(IBITPIX.EQ.16) GO TO 320 c TYPE *, IBITPIX, ' WRONG, MUST BE 16 BIT DATA' GO TO 950 C COMPRESSION FLAG. 110 IF(IDATA.EQ.1) GO TO 320 c TYPE *, 'COMPRESSION NOT ALLOWED' GO TO 950 C COLOR FLAG. 120 IF(IDATA.LE.1) GO TO 320 c TYPE *, 'COLOR NOT ALLOWED' GO TO 950 125 continue call getifel(idata,i1dar,icount,iswbf,1) c do id = 1, icount nclef = icount 126 continue mlin = ' '// 1 ' ' call charln(i1dar,icount,nclef,mlin,ncml) mlin0 = mlin c type 128, mlin(1:ncml) c write(12,*),mlin(1:ncml) ipos = index(mlin,'Year') if (ipos.eq.1) read(mlin(ipos+5:ncml),'(i5)') iyr ipos = index(mlin,'Month') if (ipos.eq.1) read(mlin(ipos+6:ncml),'(i2)') imo ipos = index(mlin,'Day') if (ipos.eq.1) read(mlin(ipos+4:ncml),'(i2)') ida ipos = index(mlin,'time') if (ipos.eq.1) then ipos2 = index(mlin,'NIL') if (ipos2 .ne.0) then tav = 19.0 else ipos2 = index(mlin,':') if ((ipos2-ipos).eq.7) then read(mlin(ipos+5:ipos2-1),'(i2)') ihr else read(mlin(ipos+5:ipos2-1),'(i1)') ihr endif mtxt = mlin(ipos2+1:ncml) ipos = index(mtxt,':') if(ipos.eq.3) then read(mtxt(1:2),'(i2)') imn else read(mtxt(1:1),'(i1)') imn endif tav = float(ihr+8)+float(imn)/60.0 endif endif do nkey = 1, num_input_tif_items mtifiden = mtifid(nkey) idlen = len(mtifid(nkey)) ipos = index(mlin0,' ') mtifiden_in = mlin0(1:ipos) if (mtifiden_in.eq.mtifid(nkey)) then ilastchar = verify(mlin0,' ',back = .true.) if (kwd_is_string(nkey)) then mtifid0 = ' / '//mtifid(nkey) mtifid0 = adjustl(mtifid0) if ((ilastchar-ipos).gt.17) then mtifkwd_v = ''''//mlin0(ipos + 1:ilastchar)//''''//mtifid0 else istart = 20 - ilastchar + ipos - 1 istart = max0(1,istart) mtifkwd_v(1:22) = ' ' mtifkwd_v(istart:20) = ''''//mlin0(ipos + 1:ilastchar)//'''' mtifkwd_v(23:70) = mtifid0 endif else mtifkwd_v0 = mlin0(ipos + 1:ilastchar) mtifkwd_v0 = adjustr(mtifkwd_v0) mtifkwd_v = mtifkwd_v0//' / '//mtifid(nkey) endif mtifkwd_value(nkey) = mtifkwd_v endif enddo 128 format(A) if(nclef.gt.0) go to 126 c type '(i4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2,f10.2)', c 1 iyr,imo,ida,ihr,imn,tav call calend(imo,ida,iyr-1900,tav,XDAY,FRACYR) ihoward = floor(xday) tav = (xday - ihoward)*24.0 call ephem2(ihoward, tav ,1.0) sun_earth_fac = fe/4.65362E-03 c type *,ihoward,sun_earth_fac c enddo go to 320 C ADDRESS(S) OF IMAGE DATA. 130 NADI = ICOUNT IF(ITYP.EQ.3) GO TO 150 C INTEGER*4 ADDRESS. IF(IVAF.NE.0) GO TO 140 IMAD(1) = IDATA !ONE ADDRESS ONLY. GO TO 320 140 IF(NADI.GT.400) GO TO 170 C GET NADI 32 BIT ADDRESSES INTO IMAD. CALL GETIFEL(IDATA, IMAD, ICOUNT*4, ISWBF, 1) GO TO 320 C INTEGER*2 ADDRESS. 150 IF(IVAF.NE.0) GO TO 160 C 1 OR 2. IMAD(1) = IREF(1) !IREF EQUIVALENCED TO IDATA. IMAD(2) = IREF(2) GO TO 320 160 IF(NADI.LE.400) GO TO 180 C TOO MANY ADDRESSES FOR BUFFER. 170 TYPE *, 'IMAGE IN TOO MANY PIECES, NEED MORE SPACE' GO TO 950 180 CALL GETIFEL(IDATA, IMAD, ICOUNT*2, ISWBF, 1) C RE-CONSTITUTE I*2 ADDRESSES INTO I*4 ARRAY. ILV = ICOUNT*2 - 1 DO 190 J=ICOUNT,1,-1 IMAD(J) = 0 I2MAD(ILV) = I2MAD(J) !I2MAD EQUIVALENCED TO IMAD ILV = ILV - 2 190 CONTINUE GO TO 320 C NUMBER OF VALUES PER PIXEL. 200 IF(IDATA.EQ.1) GO TO 320 c TYPE *, 'ONLY 1 VALUE PER PIXEL ALLOWED, COLOR?' GO TO 950 C NO. OF ROWS PER STRIP. 210 NROWPS = IDATA GO TO 320 C X PIXELS PER UNIT (INCH). 220 sclPIX1 = RDATA GO TO 320 C Y PIXELS PER UNIT. 230 sclPIX2 = RDATA GO TO 320 C UNITS FLAG. 240 GO TO (250, 260, 270), IDATA 250 MRESU = '??? ' !UNKNOWN UNITS. GO TO 320 260 MRESU = 'INCH' GO TO 320 270 MRESU = 'CM. ' 320 CONTINUE IREBF = 0 NAX1A = NAXIS1 IF(IASPEC.NE.0) NAX1A = NAXIS1*6/5 NAX2A = NAXIS2 NELO = NAX1A*NAX2A IF(NELO.GE.4000000.AND.NELO.LE.9000000) GO TO 321 C IF ARRAY SIZE IS .GT. 3000**2 OR .LT. 2000**2, PREPARE TO REBIN TO C ABOUT 2600**2 . c FAC = DSQRT(6680000.D0/DFLOAT(NELO)) c FAC = 1.0 fac = 0.595 NAX1A = FLOAT(NAXIS1)*FAC nax1a = 3*(nax1a/3) c IF(IASPEC.NE.0) NAX1A = FLOAT(NAXIS1)*FAC*6.0/5.0 + 0.5 NAX2A = FLOAT(NAXIS2)*FAC nax2a = 3*(nax2a/3) fac = nax1a/naxis1 sclpix1 = sclpix1/fac sclpix2 = sclpix2/fac IREBF = 1 c IREBF = 0 !TURN OFF REBINNING c TYPE *, 'NEW NX,NY ', NAX1A, NAX2A, FAC C SECTION TO FORMAT FITS HEADER. 321 continue c NDXO = 0 NBTR = NAXIS1*2 INIT = 1 ILC = 1 C WRITE FITS HEADER. IREC = 1 c WRITE(3,REC=IREC) BUFO IREC = IREC + 1 c c c C REFORMAT TIF DATA AND WRITE TO FITS FILE. DO 430 L=1,NADI ILINA = IMAD(L) !TIF IMAGE LINE ADDRESS. IF (NROWPS.GT.NAXIS2) NROWPS = NAXIS2 DO 420 K=1,NROWPS CALL GETIFEL(ILINA, IBUF, NBTR, ISWBF, 1) ILINA = ILINA + NBTR IF(IREBF.EQ.0) GO TO 416 C NEED TO REBIN, FLOAT A LINE. DO 410 J=1,NAXIS1 ARI(J) = IFOUR(IBUF(J)) 410 CONTINUE 411 CALL WREBIN(ARI, NAXIS1, NAXIS2, INIT, ARO, WEI, 1 NAX1A, NAX2A, NSO) IF(NSO.LE.0) GO TO 420 DO 414 J=1,NSO DO 412 I=1,NAX1A INT = ARO(I) IF(IFLIP.EQ.0) THEN INT = (INT + 1)/2 ELSE INT = (65535 - INT + 1)/2 ENDIF IF(INT.GT.32767) INT = 32767 NDXO = NDXO + 1 IPIC(NDXO) = INT 412 CONTINUE ILC = ILC + 1 IF(MOD(ILC,500).EQ.0) THEN ANT = (FLOAT(ILC)/FLOAT(NAX2A))*100. TYPE 413, ANT 413 FORMAT(1H+,' REBINNING',F8.2,' %') ENDIF 414 CONTINUE IF(ILC.LE.NAX2A) GO TO 411 GO TO 432 C NO REBINNING NEEDED, SAVE ARRAY IN MEMORY. 416 DO 417 J=1,NAX1A IF(IFLIP.EQ.0) THEN INT = (IFOUR(IBUF(J)) + 1)/2 ELSE INT = (65535 - IFOUR(IBUF(J)) + 1)/2 ENDIF IF(INT.GT.32767) INT = 32767 NDXO = NDXO + 1 IPIC(NDXO) = INT 417 CONTINUE ILC = ILC + 1 IF (ILC.GT.NAX2A) GO TO 432 420 CONTINUE 430 CONTINUE C DO A HISTOGRAM, SAVE IN FILE TIFAC. 432 CONTINUE c The following set up the dynamic storage arrays used for the radius c determination. if (allocated) then if ((nax1a.ne.nax1asave).or.(nax2a.ne.nax2asave)) then allocated = .false. deallocate(ARR) deallocate(avearr) deallocate(difarr) deallocate(mask) deallocate(arr0) deallocate(xave) deallocate(yave) deallocate(xslopeave) deallocate(yslopeave) deallocate(slopeave) endif endif if (.not.allocated) then ALLOCATE(ARR(1:NAX1A,1:NAX2A)) ALLOCATE(avearr(1:NAX1A,1:NAX2A)) ALLOCATE(difarr(1:NAX1A,1:NAX2A)) allocate(mask(1:nax1a,1:nax2a)) allocate(arr0(1:nax1a,1:nax2a)) allocate(xave(1:nax2a)) allocate(yave(1:nax1a)) allocate(xslopeave(1:nax2a)) allocate(yslopeave(1:nax1a)) allocate(slopeave(1:nax1a/2)) allocated = .true. nax1asave = nax1a nax2asave = nax2a endif NBIN = 4096 VMIN = 0. VMAX = 32767. CALL HISTIN(ARO, NBIN, VMIN, VMAX, 1) DO 450 J=1,NAX2A ILI = (J - 1)*NAX1A DO 440 I=1,NAX1A ARI(I) = IPIC(ILI + I) ARR(I,J) = ARI(I) 440 CONTINUE CALL HISTDO(ARI, NAX1A, ARO, NBIN, 1) 450 CONTINUE EM = FLOAT(NBIN)/(VMAX - VMIN) BE = 1.0 - VMIN*EM DO 460 I=1,NBIN VAL = (FLOAT(I) - BE)/EM c WRITE(2,*) VAL, ARO(I) 460 CONTINUE c The following section has been moved below the determination of c the radius and the removal of the dust. C DUMP IPIC ARRAY ON DISK IN FITS FORMAT. c 600 K = 1 c DO 620 L=1,NDXO c IF(K.LE.1440) GO TO 610 c CALL RESWAB(IBFO, 1440) c WRITE(3,REC=IREC) IBFO c IREC = IREC + 1 c K = 1 c 610 IBFO(K) = IPIC(L) c K = K + 1 c 620 CONTINUE c IF(K.GT.1440) GO TO 640 c DO 630 I=K,1440 c IBFO(I) = NUL c 630 CONTINUE c CALL RESWAB(IBFO, 1440) c WRITE(3,REC=IREC) IBFO c 640 CLOSE(UNIT=3) c TYPE *, 'NO. REC', IREC c Beginning of the section to remove dust and pits. c mask is the array where the number of times a pixel is corrected of dust and c pits is summed. The array starts in arr(i,j) and is moved successively c to avearr(i,j) and difarr(i,j). In all these corrections a series of c offset arrays are used. These consist of offsets in i and j stored as pairs c of type pr. The index k runs through all the pairs and the set of offsets c is arbitrary being designed to carry our various tasks. The first of these c is ij(k) which gives a rough circle of radius 2 pixels. The deviation between c a pixel and the average over its surrounding circle gives a measure of the c discordance of the pixel. If the deviation exceeds a limit, the value of c array at the pixel is replaced by its average. This replacement is vetoed c if any of the pixels in the circle are more than -3 times the deviation of c the target pixel. type *, 'Removing dust and pits.' ireplace = 0 do j = 1, nax2a do i = 1, nax1a mask(i,j) = 0 arr0(i,j) = arr(i,j) end do end do c The replacement tests are done four times with a decreasing deviation limit. rep: do irep = 1, 4 select case(irep) case (1) limit= 1000 case (2) limit = 800 case (3:4) limit = 700 end select Dfj1: do j = 1, nax2a Dfi1: do i = 1, nax1a c Do not consider pixels within 2 of any edge. if ((i.le.2).or.(j.le.2).or. 1 (i.ge.nax1a-1).or.(j.ge.nax2a-1)) then avsum = arr(i,j) else avsum = 0.0 c This loop gives the average in the surrounding circle. do k = 1, nptave avsum = avsum + float(arr(i+ij(k)%ii,j+ij(k)%jj)) end do avsum = avsum/nptave end if avearr(i,j) = anint(avsum) end do dfi1 end do dfj1 c Calculate the deviations between the pixel and the corresponding circle average. Dfj2: do j = 1, nax2a Dfi2: do i = 1, nax1a difarr(i,j) = arr(i,j) - avearr(i,j) end do dfi2 end do dfj2 c Determine if the replacement should be vetoed. Dfj3: do j = 3, nax2a-2 Dfi3: do i = 3, nax1a-2 c The veto test is done separately for positive and negative deviations. if (difarr(i,j).gt.limit) then c This is the veto limit for positive deviations. itest= -3*difarr(i,j) good = .true. k =1 test1:do while (good.and.(k.le.nptave)) good = (difarr(i+ij(k)%ii,j+ij(k)%jj).gt.itest) k = k+1 end do test1 if (good) then arr(i,j) = avearr(i,j) ireplace=ireplace+1 mask(i,j) = mask(i,j) +1 endif else if (difarr(i,j).lt.-limit) then c This is the veto limit for negative deviations. itest= - 3*difarr(i,j) good = .true. k =1 test2:do while (good.and.(k.le.nptave)) good = (difarr(i+ij(k)%ii,j+ij(k)%jj).lt.itest) k = k+1 end do test2 if (good) then c The point is discordant and not near a more discordant point, do replacement. arr(i,j) = avearr(i,j) ireplace=ireplace+1 mask(i,j) = mask(i,j) + 1 endif endif end do dfi3 end do dfj3 type *, 'No. points replaced:', ireplace write(8,*) 'Num. dusty points replaced:' , ireplace end do rep c This is the end of the dust and pits iterative replacement loop. c Next calculate the deviation relative to a circle having a radius c of 4 pixels along with the rms variation of this deviation over c a 5 by 5 array of points centered on the target point. c This deviation is an indication of image variation due to solar features c being present on the image. It is only used to find the first guess of the c center of the solar image. type *, 'Getting rms image for first guess center position.' Dfj4: do j = 1, nax2a Dfi4: do i = 1, nax1a if ((i.le.4).or.(j.le.4).or. 1 (i.ge.nax1a-4).or.(j.ge.nax2a-4)) then avsum = arr(i,j) else avsum = 0.0 do k = 1, nptrad c This is the sum over the radius 4 pixel circle. avsum = avsum + float(arr(i+ijr(k)%ii,j+ijr(k)%jj)) end do avsum = avsum/nptrad end if avearr(i,j) = anint(avsum) end do dfi4 end do dfj4 Dfj5: do j = 1, nax2a Dfi5: do i = 1, nax1a C This array has the deviation from the 4 pixel circular average. difarr(i,j) = -(arr(i,j) - avearr(i,j)) end do dfi5 end do dfj5 Dfj6: do j = 1, nax2a Dfi6: do i = 1, nax1a if ((i.le.2).or.(j.le.2).or. 1 (i.ge.nax1a-2).or.(j.ge.nax2a-2)) then avsum = 0.0 else avsum = 0.0 do k = 1, nptms c Here is the rms variation due to solar features. avsum = avsum + float(difarr(i+ijms(k)%ii,j+ijms(k)%jj)**2) end do avsum = avsum/nptms end if avearr(i,j) = anint(sqrt(avsum)) end do dfi6 end do dfj6 c Find an estimate for the rms variation. nptsum = 0.0 sumwt = 0.0 sumval = 0.0 do i = nax1a/4, 3*(nax1a/4) do j = nax2a/4, 3*(nax2a/4) sumwt = sumwt + avearr(i,j) sumval = sumval + arr(i,j) nptsum = nptsum + 1 enddo enddo amsave = sumwt/nptsum valave = sumval/nptsum write(8,*) 'Average ms: ',amsave c These next two loops get the first guess for the image center based on the c assumption that the plate background and scattered light are smooth. Use a c weighted average pixel sum to estimate the image center. sumsumx = 0.0 nptxsum = 0 avx: do j = 1, nax2a sumwtx = 0.0 sumwt = 0.0 avx1:do i = 1, nax1a if (avearr(i,j).lt.ams_fac*amsave) then wt = 0.0 else wt = 1.0 endif sumwtx =sumwtx + wt*i sumwt = sumwt + wt end do avx1 if (sumwt.gt.0.0) then xave(j) = sumwtx/sumwt else xave(j) = nax1a/2 endif if ((j.gt.nax2a/10).and.(j.lt.(nax2a/10)*9)) then sumsumx = sumsumx + xave(j)*sumwt nptxsum = nptxsum + sumwt end if end do avx x0 = sumsumx/nptxsum ix0 = nint(x0) c The following printout was for testing purposes. c OPEN(UNIT=4, FILE='g:\dump\xave.txt', STATUS='REPLACE', c 1 CARRIAGECONTROL='LIST') c write(4,*) 'X0= ',x0 c do j = 1, nax2a c write(4,*) j,xave(j) c end do c close(unit=4) sumsumy = 0.0 nptysum = 0 avy: do i = 1, nax1a sumwty = 0.0 sumwt = 0.0 avy1:do j = 1, nax2a if (avearr(i,j).lt.ams_fac*amsave) then wt = 0.0 else wt = 1.0 endif sumwty =sumwty + wt*j sumwt = sumwt + wt end do avy1 if (sumwt.ne.0.0) then yave(i) = sumwty/sumwt else yave(i) = nax2a/2 endif if ((i.gt.nax1a/10).and.(i.lt.(nax1a/10)*9)) then sumsumy = sumsumy + yave(i)*sumwt nptysum = nptysum + sumwt end if end do avy y0 = sumsumy/nptysum jy0 = nint(y0) c Another test printout. c OPEN(UNIT=4, FILE='g:\dump\yave.txt', STATUS='REPLACE', c 1 CARRIAGECONTROL='LIST') c write(4,*) 'Y0= ',y0 c do i = 1, nax1a c write(4,*) i, yave(i) c end do c close(unit=4) type *, 'First trial center at: X0, Y0: ',X0,Y0 write(8,*) 'First trial center at: X0, Y0: ',X0,Y0 c Now find the solar edge using a smoothed image gradient. The offset table c ijrow is used for the portions of the image either near verticle or near c horizontal. In each case averages over the pixels exterior to the point c are subtracted from corresponding pixels interior to the point. For ijrow c the array is symmetrically placed relative to the target pixel. type *, 'Calculate the radial slope-gram.' slopesum = 0.0 nptslope = 0 slopemax = 0.0 slopemin = 32000.0 slopexpsum = 0.0 nptslopexp = 0 slopexpmax = 0.0 slopexmsum = 0.0 nptslopexm = 0 slopexmmax = 0.0 slopeypsum = 0.0 nptslopeyp = 0 slopeypmax = 0.0 slopeymsum = 0.0 nptslopeym = 0 slopeymmax = 0.0 Dfj7: do j = 1, nax2a Dfi7: do i = 1, nax1a dx = i - x0 dy = j - y0 idx = i - ix0 jdy = j - jy0 if ((i.le.9).or.(j.le.9).or. 1 (i.ge.nax1a-9).or.(j.ge.nax2a-9).or. c Do not consider slopes in the inner 25 percent of the sun's apparent c disk. 1 ((iabs(idx).le.nint(nax1a/4.5)) 2 .and.(iabs(jdy).le.nint(nax2a/4.5) ))) then dfsum = 0.0 reject = .true. else reject = .false. c choose vertical or horizontal: ijswp is true close to the y axis. ijswp = (iabs(jdy).ge.iabs(idx)) c use the diagonal case if the slope is near 45 degrees than dy/dx = 0.6 c or dx/dy = - 0.6 . if (ijswp) then ijdiagonal = (3*abs(dy).lt. 5*abs(dx)) else ijdiagonal = (5*abs(dy).gt.3*abs(dx)) endif c pick the quadrant from the signs of idx and idy. if (idx.gt.0) then isn = 1 else isn = -1 end if if (jdy.gt.0) then jsn = 1 else jsn = -1 endif dfsum = 0.0 c Calculate the average gradient at the point from the difference between the c averages over two blobs oriented so that c the positive contribution is closer to the image center than the negative c contribution. diag:if (ijdiagonal) then npts = nptdiag do k = 1, nptdiag dfsum = dfsum 1 + float(arr(i-isn*ijdiag(k)%ii,j-jsn*ijdiag(k)%jj)) 1 - float(arr(i+isn*ijdiag(k)%ii,j+jsn*ijdiag(k)%jj)) end do else npts = nptrow swap:if (ijswp) then do k = 1, nptrow dfsum=dfsum 1 + float(arr(i+isn*ijrow(k)%jj,j-jsn*ijrow(k)%ii)) 1 - float(arr(i+isn*ijrow(k)%jj,j+jsn*ijrow(k)%ii)) end do else do k = 1, nptrow dfsum=dfsum 1 + float(arr(i-isn*ijrow(k)%ii,j+jsn*ijrow(k)%jj)) 1 - float(arr(i+isn*ijrow(k)%ii,j+jsn*ijrow(k)%jj)) end do end if swap end if diag dfsum = dfsum/(2*npts) end if difarr(i,j) = anint(dfsum) c Add up all the slopes to compute the average of the absolute value of c the slope over the whole array excluding those parts near disk center c and those too close to the edge of the array. This average is used as c the scaling factor to define the edge of the solar image as the point c where the average over a circular arc centered on the trial disk center c exceeds a factor times the global average slope. if (.not. reject) then slopesum = slopesum + abs(dfsum) nptslope = nptslope + 1 if (difarr(i,j).gt.slopemax) slopemax = difarr(i,j) if (abs(difarr(i,j)).lt.slopemin) slopemin = difarr(i,j) if (difarr(i,j).ge.0.0) then if ((idx.gt.0).and.(iabs(jdy).lt.nax2a/10)) then slopexpsum = slopexpsum + abs(dfsum) nptslopexp = nptslopexp + 1 if (difarr(i,j).gt.slopexpmax) slopexpmax = difarr(i,j) else if ((idx.lt.0).and.(iabs(jdy).lt.nax2a/10)) then slopexmsum = slopexmsum + abs(dfsum) nptslopexm = nptslopexm + 1 if (difarr(i,j).gt.slopexmmax) slopexmmax = difarr(i,j) else if ((jdy.gt.0).and.(iabs(idx).lt.nax1a/10)) then slopeypsum = slopeypsum + abs(dfsum) nptslopeyp = nptslopeyp + 1 if (difarr(i,j).gt.slopeypmax) slopeypmax = difarr(i,j) else if ((jdy.lt.0).and.(iabs(idx).lt.nax1a/10)) then slopeymsum = slopeymsum + abs(dfsum) nptslopeym = nptslopeym + 1 if (difarr(i,j).gt.slopeymmax) slopeymmax = difarr(i,j) end if end if end if end do dfi7 end do dfj7 offmn= offset_distance_minimum ! 0.1 offmx= offset_distance_maximum ! 0.4 ioffmn=nint(offmn*rzero*sun_earth_fac) ioffmx=nint(offmx*rzero*sun_earth_fac) cmn = sqrt(1.0-offmn**2) cmx = sqrt(1.0-offmx**2) cs = line_veto !0.3 dstv = dust_veto !10.0 ef = edge_fac aveslope = slopesum/nptslope aveslopexp = slopexpsum/nptslopexp aveslopexm = slopexmsum/nptslopexm aveslopeyp = slopeypsum/nptslopeyp aveslopeym = slopeymsum/nptslopeym c type *,'Ave grad, min grad, max grad = ', aveslope,slopemin, c 1 slopemax write(8,*) 'Ave grad, min grad, max grad = ', aveslope,slopemin, 1 slopemax write(8,*) 'xp Ave slo., max slo. = ', aveslopexp,slopexpmax write(8,*) 'xm Ave slo., max slo. = ', aveslopexm,slopexmmax write(8,*) 'yp Ave slo., max slo. = ', aveslopeyp,slopeypmax write(8,*) 'ym Ave slo., max slo. = ', aveslopeym,slopeymmax c Second round of weighted average pixel sum to improve estimate of the image center. sumsumx = 0.0 nptxsum = 0 avx2: do j = 1, nax2a sumwtx = 0.0 sumwt = 0.0 avx3:do i = 1, nax1a if((avearr(i,j).lt.ams_fac*amsave) 1 .or.(difarr(i,j).lt.0)) then wt = 0.0 else wt = 1.0 endif sumwtx =sumwtx + wt*i sumwt = sumwt + wt end do avx3 if (sumwt.gt.0.0) then xave(j) = sumwtx/sumwt else xave(j) = nax1a/2 endif if ((j.gt.nax2a/10).and.(j.lt.(nax2a/10)*9)) then sumsumx = sumsumx + xave(j)*sumwt nptxsum = nptxsum + sumwt end if end do avx2 x0 = sumsumx/nptxsum ix0 = nint(x0) sumsumy = 0.0 nptysum = 0 avy2: do i = 1, nax1a sumwty = 0.0 sumwt = 0.0 avy3:do j = 1, nax2a if((avearr(i,j).lt.ams_fac*amsave) 1 .or.(difarr(i,j).lt.0)) then wt = 0.0 else wt = 1.0 endif sumwty =sumwty + wt*j sumwt = sumwt + wt end do avy3 if (sumwt.ne.0.0) then yave(i) = sumwty/sumwt else yave(i) = nax2a/2 endif if ((i.gt.nax1a/10).and.(i.lt.(nax1a/10)*9)) then sumsumy = sumsumy + yave(i)*sumwt nptysum = nptysum + sumwt end if end do avy2 y0 = sumsumy/nptysum jy0 = nint(y0) type *, 'Second trial center at: X0, Y0: ',X0,Y0 write(8,*) 'Second trial center at: X0, Y0: ',X0,Y0 c rtestsq = (995.0*sun_earth_fac)**2 c ix00 = ix0 c jy00 = jy0 c do ixycorrect = 1 , 50 c iycount = 0 c ratypm = aveslopeyp/aveslopeym c ratxpm = aveslopexp/aveslopexm c do idi = ioffmn, ioffmx c if (iabs(idi).gt.ioffmn) then c jdj = nint(sqrt(rtestsq-float(idi)**2)) c if (difarr(ix0+idi,jy0+jdj).gt.difarr(ix0-idi,jy0-jdj)*ratypm) c 1 then c iycount = iycount + 1 c else c iycount = iycount - 1 c endif c if (difarr(ix0-idi,jy0+jdj).gt.difarr(ix0+idi,jy0-jdj)*ratypm) c 1 then c iycount = iycount + 1 c else c iycount = iycount - 1 c endif c endif c enddo c if (iycount.eq.0) then c continue c else if (iycount.gt.0) then c jy0 = jy0 + 1 c else c jy0 = jy0 - 1 c endif c ixcount = 0 c do jdj = ioffmn, ioffmx c if (iabs(jdj).gt.ioffmn) then c idi = nint(sqrt(rtestsq-float(jdj)**2)) c if (difarr(ix0+idi,jy0+jdj).gt.difarr(ix0-idi,jy0-jdj)*ratxpm) c 1 then c ixcount = ixcount + 1 c else c ixcount = ixcount - 1 c endif c if (difarr(ix0+idi,jy0-jdj).gt.difarr(ix0-idi,jy0+jdj)*ratxpm) c 1 then c ixcount = ixcount + 1 c else c ixcount = ixcount - 1 c endif c endif c enddo c if (ixcount.eq.0) then c continue c else if (ixcount.gt.0) then c ix0 = ix0 + 1 c else c ix0 = ix0 - 1 c endif c end do c write(8,*) 'Alt ix0,jy0',ix0, jy0, c 1 ' Kept ix0,jy0: ',ix00,jy00 cc type *, 'adjusted ix0,jy0',ix0, jy0 c ix0 = ix00 c jy0 = jy00 do j = 10, nax2a -10 slopesum = 0.0 nptslope = 0 idj = iabs(j - jy0) ioffavemin = ioffmn ioffavemax = ioffmx + 100 islopemin = nint(995.0*sun_earth_fac) if (idj.lt.islopemin) then ioffavemin = imax0(ioffmn, 1 nint(sqrt((rzero*sun_earth_fac)**2 - float(idj)**2))) endif do i = ix0 - ioffavemax, ix0 + ioffavemax if (iabs(i).gt.ioffavemin) then slopesum = slopesum + difarr(i,j) nptslope = nptslope + 1 endif enddo xslopeave(j) = slopesum/nptslope c write(9,*) 'j,xslopeave', j, xslopeave(j) enddo do i = 10, nax1a -10 slopesum = 0.0 nptslope = 0 idi = iabs(i - ix0) joffavemin = ioffmn joffavemax = ioffmx + 100 if (idi.lt.islopemin) then joffavemin = imax0(ioffmn, 1 nint(sqrt((rzero*sun_earth_fac)**2 - float(idi)**2))) endif do j = jy0 - joffavemax, jy0 + joffavemax if (iabs(j).gt.joffavemin) then slopesum = slopesum + difarr(i,j) nptslope = nptslope + 1 endif enddo yslopeave(i) = slopesum/nptslope c write(9,*) 'i,yslopeave',i, yslopeave(i) enddo istart_radius = nint(sqrt(rzero**2-ioffmn**2)*sun_earth_fac) i = ix0-istart_radius do while ((yslopeave(i).gt.-4.0*aveslopexm).and.(i.gt.13)) i = i-1 enddo ixmst = i i = ix0+istart_radius do while ((yslopeave(i).gt.-4.0*aveslopexp).and.(i.lt.(nax1a-13))) i = i+1 enddo ixpst = i j = jy0-istart_radius do while ((xslopeave(j).gt.-4.0*aveslopeym).and.(j.gt.13)) j = j-1 enddo jymst = j j = jy0+istart_radius do while ((xslopeave(j).gt.-4.0*aveslopeyp).and.(j.lt.(nax2a-13))) j = j+1 enddo jypst = j c type *, ixmst,ixpst,jymst,jypst write(8,*) 'ixm,ixp,jym,jyp', ixmst,ixpst,jymst,jypst c close(unit = 9) c OPEN(UNIT=4, FILE='g:\dump\radxp.txt', STATUS='REPLACE', c 1 CARRIAGECONTROL='LIST') c Start all circular arcs with a radius equal to the minimum distance c between the array edge and the trial disk center. c Scan inwards calculating the average of the slope along an arc centered c on disk center. Do this for quadrants positive and negative x and positive c and negative y. radxp=30.0 radyp=0 iraddo=0 scatxpslope = 0.0 scatxmslope = 0.0 scatypslope = 0.0 scatymslope = 0.0 scfac = scattered_light_factor ! 0.7 raddo: 1 do while (((abs(radxp-radxm).gt.1.0).or. 2 (abs(radyp-radym).gt.1.0)).and.(iraddo.lt.9)) select case (iraddo) case(0,1,2) grfac = grad_factor rstart0 = rzero*1.17 rlimit0 = rzero*1.135 rmin0 = rmin*0.95 case(3) grfac = grad_factor*1.0 rstart0 = rzero*1.145 rlimit0 = rzero*1.085 rmin0 = rmin case(4) grfac = grad_factor*1.0 rstart0 = rzero*1.10 rlimit0 = rzero*1.075 case(5) grfac = grad_factor*1.0 rstart0 = rzero*1.085 rlimit0 = rzero*1.06 case(6) grfac = grad_factor*1.0 rstart0 = rzero*1.075 rlimit0 = rzero*1.05 end select idemerit = 0 c Start of x plus section. router = (ixpst-ix0 - 11)*sqrt(1.0+offmn**2) rstart = amin1(router,rstart0*sun_earth_fac) rlimit = amin1(router,rlimit0*sun_earth_fac) radxp = nax1a-ix0-15 radxp = amin1(radxp,rstart) radxpmax = radxp slopesum = 0.0 slopevmin = 0.0 amssum = 0.0 amsxpave = 0.0 ircount=1 xpslopeave = 0.0 xpslopemax = 0.0 carry_on = .true. do while (carry_on) radxpsq = radxp**2 slopevmin = 0.0 slopesum = 0.0 yslopeavemx = 0.0 amssum = 0.0 nptslope = 0 do jrad = nint(offmn*radxp), nint(offmx*radxp) rad= jrad idi = nint(sqrt(radxpsq - rad**2)) nptslope = nptslope + 2 ipt = ix0+idi jpt1 = jy0+jrad jpt2 = jy0-jrad if ((ipt.lt.nax1a).and.(ipt.ge.1).and.(jpt1.lt.nax2a) 1 .and.jpt2.ge.1) then slopesum = slopesum + difarr(ipt,jpt1) 1 + difarr(ipt,jpt2) amssum = amssum + avearr(ipt,jpt1) + avearr(ipt,jpt2) if(difarr(ipt-7,jpt1).lt.slopevmin)slopevmin=difarr(ipt-7,jpt1) if(difarr(ipt-7,jpt2).lt.slopevmin)slopevmin=difarr(ipt-7,jpt2) if(difarr(ipt-13,jpt1).lt.slopevmin)slopevmin=difarr(ipt-13,jpt1) if(difarr(ipt-13,jpt2).lt.slopevmin)slopevmin=difarr(ipt-13,jpt2) if(difarr(ipt+7,jpt1).lt.slopevmin)slopevmin=difarr(ipt+7,jpt1) if(difarr(ipt+7,jpt2).lt.slopevmin)slopevmin=difarr(ipt+7,jpt2) if(difarr(ipt+13,jpt1).lt.slopevmin)slopevmin=difarr(ipt+13,jpt1) if(difarr(ipt+13,jpt2).lt.slopevmin)slopevmin=difarr(ipt+13,jpt2) if(yslopeave(ipt).gt.yslopeavemx) yslopeavemx = yslopeave(ipt) else c type *, 'Bad point xp: ', ipt,jpt1,jpt2 stop endif end do xpslopeave = slopesum/nptslope amsxpave = amssum/nptslope slopeave(ircount)=xpslopeave radxp = radxp - 1 ircount = ircount+1 ircountxp = ircount irc1 = max0(ircountxp - 25, 2) nirc = 0 slopesum = 0.0 do irc = irc1, ircountxp - 10 if (slopeave(irc).gt.-grfac*aveslope) then slopesum = slopesum + slopeave(irc) nirc = nirc +1 endif enddo if (nirc.gt.0) then scatxpslope0 = slopesum/nirc endif dust_test0=dust_test line_test0=line_test edge_test0=edge_test grad_test0=grad_test blank_test0=blank_test rmin_test0=rmin_test rlim_test0=rlim_test grad_test = ((xpslopeave-scfac*scatxpslope).lt.grfac*aveslopexp) blank_test = ((avearr(ix0+nint(radxp),jy0).eq.0.0) 3 .or.(avearr(ix0+nint(radxp)+13,jy0).eq.0.0)) ys1 = yslopeave(ix0+nint(cmn*radxp)-7) ys2 = yslopeave(ix0+nint(cmn*radxp)-13) ys3 = yslopeave(ix0+nint(cmn*radxp)+7) ys4 = yslopeave(ix0+nint(cmn*radxp)+13) line_test = 1 ((yslopeave(ix0+nint(cmn*radxp)-7).lt.-cs*grfac*aveslopexp) 2.or.(yslopeave(ix0+nint(cmn*radxp)-13).lt. 3 -cs*grfac*aveslopexp) 2.or.(yslopeave(ix0+nint(cmn*radxp)+7).lt.-cs*grfac*aveslopexp) 2.or.(yslopeave(ix0+nint(cmn*radxp)+13).lt. 3 -cs*grfac*aveslopexp)) rmin_test = (radxp.gt.rmin0*rzero*sun_earth_fac) rlim_test = (radxp.gt.rlimit) dust_test = (slopevmin.lt.-dstv*aveslopexp) edge_test = .false. !(yslopeavemx.gt.ef*aveslopexp) carry_on = ((grad_test.or.rlim_test 1 .or.(iraddo.le.1).or.line_test.or.dust_test.or.blank_test 2 .or.edge_test).and.rmin_test) if ((.not.(dust_test.or.line_test)) 1 .and.(xpslopeave.gt.xpslopemax) 2 .and.(iraddo.le.1) 2 ) then xpslopemax = xpslopeave radxpmax = radxp endif end do if (iraddo.le.1) radxp = radxpmax xps = ' ' if (grad_test0) write(8,*) 'Grad xp' if (line_test0) write(8,*) 'Line xp' if (dust_test0) write(8,*) 'Dust xp' if (edge_test0) write(8,*) 'Edge xp' if (blank_test0) write(8,*) 'Blank xp' if (grad_test0) xps = ' G' if (line_test0) then xps = ' L' idemerit = idemerit + 1 endif if (edge_test0) then xps = ' E' idemerit = idemerit + 1 endif if (grad_test0.and.line_test0) xps = ' l' if (grad_test0.and.edge_test0) xps = ' e' if (dust_test0) then xps = ' D' idemerit = idemerit + 1 endif if (blank_test0) then xps = ' B' idemerit = idemerit + 1 endif if (rlim_test0) then xps = ' R' idemerit = idemerit + 5 endif if (.not.rmin_test) then xps = ' r' idemerit = idemerit + 5 endif if (aveslopexp.lt.(scmax*xpslopemax)) 1 aveslopexp = scmax*xpslopemax scatxpslope = amin1(scatxpslope0,0.5*xpslopeave) if (scatxpslope.lt.scatxpslope0) then xps(1:1) = 'S' idemerit = idemerit + 2 endif write(8,*) 'xp',radxp,xpslopeave,amsxpave,slopevmin c type *, 'Radius, X plus: ',radxp c Start of x minus section radxm = ix0-15 router = (ix0-ixmst - 11)*sqrt(1.0+offmn**2) rstart = amin1(router,rstart0*sun_earth_fac) rlimit = amin1(router,rlimit0*sun_earth_fac) radxm = amin1(radxm,rstart) radxmmax = radxm slopesum = 0.0 amssum = 0.0 amsxmave = 0.0 slopevmin = 0.0 xmslopeave = 0.0 xmslopemax = 0.0 carry_on = .true. ircount = 1 do while (carry_on) slopesum= 0.0 nptslope = 0 amssum = 0.0 yslopeavemx = 0.0 radxmsq = radxm**2 slopevmin =0.0 do jrad = nint(offmn*radxm), nint(offmx*radxm) rad=jrad idi = nint(sqrt(radxmsq - rad**2)) nptslope = nptslope + 2 ipt = ix0-idi jpt1 = jy0+jrad jpt2 = jy0-jrad if ((ipt.lt.nax1a).and.(ipt.ge.1).and.(jpt1.lt.nax2a) 1 .and.jpt2.ge.1) then slopesum = slopesum + difarr(ipt,jpt1) 1 + difarr(ipt,jpt2) amssum = amssum + avearr(ipt,jpt1) + avearr(ipt,jpt2) if(difarr(ipt-7,jpt1).lt.slopevmin)slopevmin=difarr(ipt-7,jpt1) if(difarr(ipt-7,jpt2).lt.slopevmin)slopevmin=difarr(ipt-7,jpt2) if(difarr(ipt-13,jpt1).lt.slopevmin)slopevmin=difarr(ipt-13,jpt1) if(difarr(ipt-13,jpt2).lt.slopevmin)slopevmin=difarr(ipt-13,jpt2) if(difarr(ipt+7,jpt1).lt.slopevmin)slopevmin=difarr(ipt+7,jpt1) if(difarr(ipt+7,jpt2).lt.slopevmin)slopevmin=difarr(ipt+7,jpt2) if(difarr(ipt+13,jpt1).lt.slopevmin)slopevmin=difarr(ipt+13,jpt1) if(difarr(ipt+13,jpt2).lt.slopevmin)slopevmin=difarr(ipt+13,jpt2) if(yslopeave(ipt).gt.yslopeavemx) yslopeavemx = yslopeave(ipt) else c type *, 'Bad point xm: ', ipt,jpt1,jpt2 stop c write(4,*) ipt,jpt1,jpt2 endif end do xmslopeave = slopesum/nptslope amsxmave = amssum/nptslope slopeave(ircount) = xmslopeave radxm = radxm - 1 ircount = ircount + 1 ircountxm = ircount irc1 = max0(ircountxm - 25, 2) nirc = 0 slopesum = 0.0 do irc = irc1, ircountxm - 10 if (slopeave(irc).gt.-grfac*aveslope) then slopesum = slopesum + slopeave(irc) nirc = nirc +1 endif enddo if (nirc.gt.0) then scatxmslope0 = slopesum/nirc endif dust_test0=dust_test line_test0=line_test edge_test0=edge_test grad_test0=grad_test blank_test0=blank_test rmin_test0=rmin_test rlim_test0=rlim_test grad_test = ((xmslopeave-scfac*scatxmslope).lt.grfac*aveslopexm) rmin_test = (radxm.gt.rmin0*rzero*sun_earth_fac) rlim_test = (radxm.gt.rlimit) dust_test = (slopevmin.lt.-dstv*aveslopexm) blank_test = ((avearr(ix0-nint(radxm),jy0).eq.0.0) 3 .or.(avearr(ix0-nint(radxm)-13,jy0).eq.0.0)) line_test = 1 ((yslopeave(ix0-nint(cmn*radxm)-7).lt.-cs*grfac*aveslopexm) 2.or.(yslopeave(ix0-nint(cmn*radxm)-13).lt. 3 -cs*grfac*aveslopexm) 2.or.(yslopeave(ix0-nint(cmn*radxm)+7).lt.-cs*grfac*aveslopexm) 2.or.(yslopeave(ix0-nint(cmn*radxm)+13).lt. 3 -cs*grfac*aveslopexm)) edge_test = .false. !(yslopeavemx.gt.ef*aveslopexm) carry_on = ((grad_test.or.rlim_test 1 .or.(iraddo.le.1).or.line_test.or.dust_test.or.blank_test 2 .or.edge_test).and.rmin_test) if ((.not.(dust_test.or.line_test)) 1 .and.(xmslopeave.gt.xmslopemax) 2 .and.(iraddo.le.1) 2 ) then xmslopemax = xmslopeave radxmmax = radxm endif end do if (iraddo.le.1) radxm = radxmmax if (grad_test0) write(8,*) 'Grad xm' if (line_test0) write(8,*) 'Line xm' if (dust_test0) write(8,*) 'Dust xm' if (edge_test0) write(8,*) 'Edge xm' if (blank_test0) write(8,*) 'Blank xm' xms = ' ' if (grad_test0) xms = ' G' if (line_test0) then xms = ' L' idemerit = idemerit + 1 endif if (edge_test0) then xms = ' E' idemerit = idemerit + 1 endif if (grad_test0.and.line_test0) xms = ' l' if (grad_test0.and.edge_test0) xms = ' e' if (dust_test0) then xms = ' D' idemerit = idemerit + 1 endif if (blank_test0) then xms = ' B' idemerit = idemerit + 1 endif if (rlim_test0) then xms = ' R' idemerit = idemerit + 5 endif if (.not.rmin_test) then xms = ' r' idemerit = idemerit + 5 endif if (aveslopexm.lt.(scmax*xmslopemax)) 1 aveslopexm = scmax*xmslopemax scatxmslope = amin1(scatxmslope0,0.5*xmslopeave) if (scatxmslope.lt.scatxmslope0) then xms(1:1) = 'S' idemerit = idemerit + 2 endif write(8,*) 'xm',radxm,xmslopeave,amsxmave,slopevmin c type *, 'Radius, X minus: ',radxm c Start of y plus section radyp = nax2a-jy0-15 router = (jypst-jy0 - 11)*sqrt(1.0+offmn**2) rstart = amin1(router,rstart0*sun_earth_fac) rlimit = amin1(router,rlimit0*sun_earth_fac) radyp = amin1(radyp,rstart) radypmax = radyp slopesum = 0.0 nptslope = 1 slopehmin = 0.0 amssum = 0.0 amsypave = 0.0 ypslopeave = 0.0 ypslopemax = 0.0 carry_on = .true. ircount = 1 do while (carry_on) slopesum=0.0 xslopeavemx = 0.0 amssum = 0.0 nptslope = 0 slopehmin =0.0 radypsq = radyp**2 do irad = nint(offmn*radyp), nint(offmx*radyp) rad= irad jdj = nint(sqrt(radypsq - rad**2)) nptslope = nptslope + 2 jpt = jy0+jdj ipt1 = ix0+irad ipt2 = ix0-irad if ((jpt.lt.nax2a).and.(jpt.ge.1).and.(ipt1.lt.nax1a) 1 .and.ipt2.ge.1) then slopesum = slopesum + difarr(ipt1,jpt) 1 + difarr(ipt2,jpt) amssum = amssum + avearr(ipt1,jpt) + avearr(ipt2,jpt) if(difarr(ipt1,jpt-7).lt.slopehmin)slopehmin=difarr(ipt1,jpt-7) if(difarr(ipt2,jpt-7).lt.slopehmin)slopehmin=difarr(ipt2,jpt-7) if(difarr(ipt1,jpt-10).lt.slopehmin)slopehmin=difarr(ipt1,jpt-10) if(difarr(ipt2,jpt-10).lt.slopehmin)slopehmin=difarr(ipt2,jpt-10) if(difarr(ipt1,jpt-13).lt.slopehmin)slopehmin=difarr(ipt1,jpt-13) if(difarr(ipt2,jpt-13).lt.slopehmin)slopehmin=difarr(ipt2,jpt-13) if(difarr(ipt1,jpt+7).lt.slopehmin)slopehmin=difarr(ipt1,jpt+7) if(difarr(ipt2,jpt+7).lt.slopehmin)slopehmin=difarr(ipt2,jpt+7) if(difarr(ipt1,jpt+13).lt.slopehmin)slopehmin=difarr(ipt1,jpt+13) if(difarr(ipt2,jpt+13).lt.slopehmin)slopehmin=difarr(ipt2,jpt+13) if(xslopeave(jpt).gt.xslopeavemx) xslopeavemx = xslopeave(jpt) else c type *, 'Bad point yp: ', ipt1,ipt2,jpt stop c write(4,*) ipt1,ipt2,jpt endif end do ypslopeave = slopesum/nptslope amsypave = amssum/nptslope slopeave(ircount) = ypslopeave radyp = radyp - 1 ircount = ircount + 1 ircountyp = ircount irc1 = max0(ircountyp - 25, 2) nirc = 0 slopesum = 0.0 do irc = irc1, ircountyp - 10 if (slopeave(irc).gt.-grfac*aveslope) then slopesum = slopesum + slopeave(irc) nirc = nirc +1 endif enddo if (nirc.gt.0) then scatypslope0 = slopesum/nirc endif dust_test0=dust_test line_test0=line_test edge_test0=edge_test grad_test0=grad_test blank_test0=blank_test rmin_test0=rmin_test rlim_test0=rlim_test grad_test =((ypslopeave-scfac*scatypslope).lt.grfac*aveslopeyp) line_test = 1 ((xslopeave(jy0+nint(cmn*radyp)-7).lt.-cs*grfac*aveslopeyp) 2.or.(xslopeave(jy0+nint(cmn*radyp)-13).lt. 3 -cs*grfac*aveslopeyp) 2.or.(xslopeave(jy0+nint(cmn*radyp)+7).lt.-cs*grfac*aveslopeyp) 2.or.(xslopeave(jy0+nint(cmn*radyp)+13).lt. 3 -cs*grfac*aveslopeyp) 2.or.(xslopeave(jy0+nint(cmn*radyp)+4).lt. 3 -cs*grfac*aveslopeyp)) blank_test = ((avearr(ix0,jy0+nint(radyp)).eq.0.0) 3 .or.(avearr(ix0,jy0+nint(radyp)+13).eq.0.0)) rmin_test = (radyp.gt.rmin0*rzero*sun_earth_fac) rlim_test = (radyp.gt.rlimit) dust_test = (slopehmin.lt.-dstv*aveslopeyp) edge_test = ((xslopeavemx.gt.ef*aveslopeyp).and.(irebf.eq.1)) carry_on = ((grad_test.or.rlim_test 1 .or.(iraddo.le.1).or.line_test.or.dust_test.or.blank_test 2 .or.edge_test).and.rmin_test) if ((.not.(dust_test.or.line_test)) 1 .and.(ypslopeave.gt.ypslopemax) 2 .and.(iraddo.le.1) 2 ) then ypslopemax = ypslopeave radypmax = radyp endif end do if (iraddo.le.1) radyp = radypmax xsa1=xslopeave(jy0+nint(cmx*radyp)-7) xsa2=xslopeave(jy0+nint(cmx*radyp)-13) xsa3=xslopeave(jy0+nint(cmx*radyp)+7) xsa4=xslopeave(jy0+nint(cmx*radyp)+13) if (grad_test0) write(8,*) 'Grad yp' if (line_test0) write(8,*) 'Line yp' if (edge_test0) write(8,*) 'Edge yp' if (dust_test0) write(8,*) 'Dust yp' if (blank_test0) write(8,*) 'Blank yp' yps = ' ' if (grad_test0) yps = ' G' if (line_test0) then yps = ' L' idemerit = idemerit + 1 endif if (edge_test0) then yps = ' E' idemerit = idemerit + 1 endif if (grad_test0.and.line_test0) yps = ' l' if (grad_test0.and.edge_test0) yps = ' e' if (dust_test0) then yps = ' D' idemerit = idemerit + 1 endif if (blank_test0) then yps = ' B' idemerit = idemerit + 1 endif if (rlim_test0) then yps = ' R' idemerit = idemerit + 5 endif if (.not.rmin_test) then yps = ' r' idemerit = idemerit + 5 endif if (aveslopeyp.lt.(scmax*ypslopemax)) 1 aveslopeyp = scmax*ypslopemax scatypslope = amin1(scatypslope0,0.5*ypslopeave) if (scatypslope.lt.scatypslope0) then yps(1:1) = 'S' idemerit = idemerit + 2 endif write(8,*) 'yp',radyp,ypslopeave,amsypave,slopehmin c type *, 'Radius, Y plus: ',radyp c Start of the y minus section radym = jy0-15 router = (jy0-jymst - 11)*sqrt(1.0+offmn**2) rstart = amin1(router,rstart0*sun_earth_fac) rlimit = amin1(router,rlimit0*sun_earth_fac) radym = amin1(radym,rstart) radymmax = radym slopesum = 0.0 amsymave = 0.0 slopehmin = 0.0 ymslopeave = 0.0 ymslopemax = 0.0 carry_on = .true. ircount = 1 do while (carry_on) slopesum= 0.0 amssum = 0.0 xslopeavemx = 0.0 nptslope = 0 radymsq = radym**2 slopehmin = 0.0 do irad = nint(offmn*radym), nint(offmx*radym) rad=irad jdj = nint(sqrt(radymsq - rad**2)) nptslope = nptslope + 2 jpt = jy0-jdj ipt1 = ix0+irad ipt2 = ix0-irad if ((jpt.lt.nax2a).and.(jpt.ge.1).and.(ipt1.lt.nax1a) 1 .and.(ipt2.ge.1)) then slopesum = slopesum + difarr(ipt1,jpt) 1 + difarr(ipt2,jpt) amssum = amssum + avearr(ipt1,jpt) + avearr(ipt2,jpt) if(difarr(ipt1,jpt-7).lt.slopehmin)slopehmin=difarr(ipt1,jpt-7) if(difarr(ipt2,jpt-7).lt.slopehmin)slopehmin=difarr(ipt2,jpt-7) if(difarr(ipt1,jpt-13).lt.slopehmin)slopehmin=difarr(ipt1,jpt-13) if(difarr(ipt2,jpt-13).lt.slopehmin)slopehmin=difarr(ipt2,jpt-13) if(difarr(ipt1,jpt+7).lt.slopehmin)slopehmin=difarr(ipt1,jpt+7) if(difarr(ipt2,jpt+7).lt.slopehmin)slopehmin=difarr(ipt2,jpt+7) if(difarr(ipt1,jpt+10).lt.slopehmin)slopehmin=difarr(ipt1,jpt+10) if(difarr(ipt2,jpt+10).lt.slopehmin)slopehmin=difarr(ipt2,jpt+10) if(difarr(ipt1,jpt+13).lt.slopehmin)slopehmin=difarr(ipt1,jpt+13) if(difarr(ipt2,jpt+13).lt.slopehmin)slopehmin=difarr(ipt2,jpt+13) if(xslopeave(jpt).gt.xslopeavemx) xslopeavemx = xslopeave(jpt) else c type *, 'Bad point ym: ', ipt1,ipt2,jpt stop c write(4,*) ipt1,ipt2,jpt endif end do ymslopeave = slopesum/nptslope amsymave = amssum/nptslope slopeave(ircount) = ymslopeave radym = radym - 1 ircount = ircount + 1 ircountym = ircount irc1 = max0(ircountym - 25, 2) nirc = 0 slopesum = 0.0 do irc = irc1, ircountym - 10 if (slopeave(irc).gt.-grfac*aveslope) then slopesum = slopesum + slopeave(irc) nirc = nirc +1 endif enddo if (nirc.gt.0) then scatymslope0 = slopesum/nirc endif dust_test0=dust_test line_test0=line_test edge_test0=edge_test grad_test0=grad_test blank_test0=blank_test rmin_test0=rmin_test rlim_test0=rlim_test grad_test = ((ymslopeave-scfac*scatymslope).lt.grfac*aveslopeym) line_test = 1 ((xslopeave(jy0-nint(cmn*radym)+7).lt.-cs*grfac*aveslopeym) 2.or.(xslopeave(jy0-nint(cmn*radym)+13).lt. 3 -cs*grfac*aveslopeym) 2.or.(xslopeave(jy0-nint(cmn*radym)-7).lt.-cs*grfac*aveslopeym) 2.or.(xslopeave(jy0-nint(cmn*radym)-13).lt. 3 -cs*grfac*aveslopeym)) blank_test = ((avearr(ix0,jy0-nint(radym)).eq.0.0) 3 .or.(avearr(ix0,jy0-nint(radym)-13).eq.0.0)) rmin_test = (radym.gt.rmin0*rzero*sun_earth_fac) rlim_test = (radym.gt.rlimit) dust_test = (slopehmin.lt.-dstv*aveslopeym) edge_test = ((xslopeavemx.gt.ef*aveslopeym).and.(irebf.eq.1)) carry_on = ((grad_test.or.rlim_test 1 .or.(iraddo.le.1).or.line_test.or.dust_test.or.blank_test 2 .or.edge_test).and.rmin_test) if ((.not.(dust_test.or.line_test)) 1 .and.(ymslopeave.gt.ymslopemax) 2 .and.(iraddo.le.1) 2 ) then ymslopemax = ymslopeave radymmax = radym endif end do if (iraddo.le.1) radym=radymmax if (grad_test0) write(8,*) 'Grad ym' if (line_test0) write(8,*) 'Line ym' if (dust_test0) write(8,*) 'Dust ym' if (edge_test0) write(8,*) 'Edge ym' if (blank_test0) write(8,*) 'Blank ym' yms = ' ' if (grad_test0) yms = ' G' if (line_test0) then yms = ' L' idemerit = idemerit + 1 endif if (edge_test0) then yms = ' E' idemerit = idemerit + 1 endif if (grad_test0.and.line_test0) yms = ' l' if (grad_test0.and.edge_test0) yms = ' e' if (dust_test0) then yms = ' D' idemerit = idemerit + 1 endif if (blank_test0) then yms = ' B' idemerit = idemerit + 1 endif if (rlim_test0) then yms = ' R' idemerit = idemerit + 5 endif if (.not.rmin_test) then yms = ' r' idemerit = idemerit + 5 endif if (aveslopeym.lt.(scmax*ymslopemax)) 1 aveslopeym = scmax*ymslopemax scatymslope = amin1(scatymslope0,0.5*ymslopeave) if (scatymslope.lt.scatymslope0) then yms(1:1) = 'S' idemerit = idemerit + 2 endif write(8,*) 'ym',radym,ymslopeave,amsymave,slopehmin c adjust image center position for the difference between the positive and c negative x and y radii. c type *, 'Radius, Y minus: ',radym if (iraddo.le.1) then write(8,*) 'ave(xp,xm,yp,ym)', 1 aveslopexp,aveslopexm,aveslopeyp,aveslopeym radxp = 100.0 + radxm ! cause convergence test to fail. else write(8,*) 'scat(xp,xm,yp,ym)', 1 scatxpslope,scatxmslope,scatypslope,scatymslope write(8,*) 'rxp,rxm,ryp,rym: ',radxp,radxm,radyp,radym x0 = x0+0.5*(radxp - radxm) y0 = y0+0.5*(radyp - radym) ix0 = nint(x0) jy0 = nint(y0) r0 = 0.25*(radxp + radxm + radyp + radym) idemerit_veto = idemerit idemerit = idemerit 1 + iabs(nint((r0 - 1020.0*sun_earth_fac)/10.0)) idemerit_rad = idemerit deltarad = (radxp + radxm - radyp - radym)/2.0 idemerit = idemerit + iabs(nint((deltarad+15.0)/10.0)) write(8,*) 'demerits, veto, rad, all: ' 1 ,idemerit_veto,idemerit_rad,idemerit xfactor = 0.5 * deltarad/r0 xstretch = 1.0 + xfactor ystretch = 1.0 - xfactor type *, 'x0, y0, r0, idm', x0, y0, r0, idemerit c type *, 'rx-ry, xstretch, ystretch = ',deltarad,xstretch,ystretch write(8,*) 'x0, y0, r0', x0, y0, r0 write(8,*) 'rx-ry, xstr, ystr = ',deltarad,xstretch,ystretch write(8,*) mtif(1:nc-1),x0,y0,r0 endif iraddo =iraddo + 1 end do raddo close(unit=8) c close(unit=9) c The section below will put into the array a circular marker indicating the determined c solar image boundary. r0sq = r0**2 ry0sq = r0sq*ystretch**2 rx0sq = r0sq*xstretch**2 do jrad = 1,nint(ystretch*r0/sqrt(2.0)) rad = jrad radsq = rad**2/ystretch**2 ijoff = nint(xstretch*sqrt(r0sq - radsq)) difarr(ix0+ijoff,jy0+jrad)= 1 mod(difarr(ix0+ijoff,jy0+jrad)+25000,28000) difarr(ix0+ijoff,jy0-jrad)= 1 mod(difarr(ix0+ijoff,jy0-jrad)+25000,28000) difarr(ix0-ijoff,jy0+jrad)= 1 mod(difarr(ix0-ijoff,jy0+jrad)+25000,28000) difarr(ix0-ijoff,jy0-jrad)= 1 mod(difarr(ix0-ijoff,jy0-jrad)+25000,28000) enddo do irad = 1,nint(xstretch*r0/sqrt(2.0)) rad = irad radsq = rad**2/xstretch**2 ijoff = nint(ystretch*sqrt(r0sq - radsq)) difarr(ix0+irad,jy0+ijoff)= 1 mod(difarr(ix0+irad,jy0+ijoff)+25000,28000) difarr(ix0-irad,jy0+ijoff)= 1 mod(difarr(ix0-irad,jy0+ijoff)+25000,28000) difarr(ix0+irad,jy0-ijoff)= 1 mod(difarr(ix0+irad,jy0-ijoff)+25000,28000) difarr(ix0-irad,jy0-ijoff)= 1 mod(difarr(ix0-irad,jy0-ijoff)+25000,28000) end do write(11,480) caseid,x0,y0,r0,deltarad,amsave,ireplace, 1 nint(aveslopexp),nint(aveslopexm), 2 nint(aveslopeyp),nint(aveslopeym), 3 nint(scatxpslope),nint(scatxmslope), 4 nint(scatypslope),nint(scatymslope),xps,xms,yps,yms,idemerit 480 format(1xa40,3f8.2,f7.1,f6.1,i7,4I5,' ',4I4,4A2,I3) c The following section uses the array visualizer to examine the array: arr c type *, 'Average Radius: ',r0,', X-Radius - Y-Radius: ', deltarad c CALL faglStartWatch (arr, ISTAT) c IF (ISTAT.NE.0) GO TO 800 c CALL faglName (arr, MTIF , ISTAT) c IF (ISTAT.NE.0) GO TO 800 c CALL faglShow (arr, ISTAT) c IF (ISTAT.NE.0) GO TO 800 c call faglEndWatch(arr, ISTAT) crpix1 = x0 crpix2 = y0 cdelt1 = 2.0/(radxp+radxm) cdelt2 = 2.0/(radyp+radym) do 745 ioutloc = 1 , 2 select case(ioutloc) case(1) nb = 1 case(2) nb = 3 end select DO 522 I=1,720 BUFO(I) = ' ' 522 CONTINUE IMBA = 1 WRITE(MLIN,530) 'SIMPLE ', 'T', 'FITS STANDARD' 530 FORMAT(A,'= ',A20,' / ',A) CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 1 IMBA = IMBA + 20 WRITE(MLIN,540) 'BITPIX ', IBITPIX, 'FITS BITS/PIXEL' 540 FORMAT(A,'= ',I20,' / ',A) CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 2 IMBA = IMBA + 20 WRITE(MLIN,540) 'NAXIS ', 2, 'NUMBER OF AXES' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 3 IMBA = IMBA + 20 WRITE(MLIN,540) 'NAXIS1 ', NAX1A/nb, 'X AXIS NO. PIX' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 4 IMBA = IMBA + 20 WRITE(MLIN,540) 'NAXIS2 ', NAX2A/nb, 'Y AXIS NO. PIX' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 5 IMBA = IMBA + 20 WRITE(MLIN,540) 'BLANK ', -32768, 'ARRAY NULL VALUE' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 6 IMBA = IMBA + 20 write(mlin,545) iyearrun,imonthrun,idayrun, 1 ihourrun,iminrun,isecrun 545 format( 1 'HISTORY Created by tiftofits version 1 on ',i4,'-' 2 ,i2.2,'-',i2.2,' at ',i2.2,':',i2.2,':',i2.2) CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 7 IMBA = IMBA + 20 write(mlin,546) 546 format( 1 'HISTORY 1-3 pixel sized dust and emulsion pits removed by Laplac 2ian filter') CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 7 IMBA = IMBA + 20 WRITE(MLIN,550) 'BSCALE ', 1.0, 'REAL = DATA*BSCALE + BZERO' 550 FORMAT(A,'= ',F20.7,' / ',A) CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 7 IMBA = IMBA + 20 WRITE(MLIN,550) 'BZERO ', 0.0, ' ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 8 IMBA = IMBA + 20 WRITE(MLIN,530) 'BUNIT ', '''INTENSITY'' ', 1 'ARBITRARY UNITS' imba_quant = imba CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 9 IMBA = IMBA + 20 WRITE(MLIN,530) 'ORIGIN ', '''MT. WILSON'' ', ' ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 10 IMBA = IMBA + 20 WRITE(MLIN,530) 'TELESCOP', '''60 FT'' ', ' ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 11 IMBA = IMBA + 20 imba_obj = imba WRITE(MLIN,530) 'OBJECT ', '''SUN'' ', ' ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 12 IMBA = IMBA + 20 IF(WAVEL.EQ.0.0) THEN WRITE(MLIN,530) 'DATA-TYP', '''DIRECT'' ', ' ' ELSE WRITE(MLIN,530) 'DATA-TYP', '''CA SPECTROHELIOGRM''', ' ' ENDIF CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 13 IMBA = IMBA + 20 WRITE(MLIN,560) 'DATE-OBS', IYR, IMO, IDA, 'YEAR-MO-DA' 560 FORMAT(A,'= ','''',I4,'-',I2.2,'-',I2.2,'''',8X,' / ',A) CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 14 IMBA = IMBA + 20 WRITE(MLIN,550) 'TAVG-UT ', TAV, 'AVERAGE TIME, U.T.' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 15 IMBA = IMBA + 20 IF(WAVEL.NE.0.0) THEN WRITE(MLIN,550) 'LAMBDA ', WAVEL, 'WAVELENGTH (ANGSTROMS)' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 16 IMBA = IMBA + 20 ENDIF WRITE(MLIN,570) 'SCALPIX1', SclPIX1/nb, 'X', MRESU 570 FORMAT(A,'= ',F20.7,' / ',A,' PIXELS PER ',A) CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 17 IMBA = IMBA + 20 WRITE(MLIN,570) 'SCALPIX2', SclPIX2/nb, 'Y', MRESU CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 18 IMBA = IMBA + 20 WRITE(MLIN,575) 'CRPIX1 ', CRPIX1/nb 575 FORMAT(A,'= ',F20.7,' / Ref point1.') CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 19 IMBA = IMBA + 20 WRITE(MLIN,576) 'CRPIX2 ', CRPIX2/nb 576 FORMAT(A,'= ',F20.7,' / Ref point2.') CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 20 IMBA = IMBA + 20 WRITE(MLIN,577) 'CDELT1 ', CDelt1*nb 577 FORMAT(A,'= ',F20.8,' / Solar Scale Axis 1.') CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 21 IMBA = IMBA + 20 WRITE(MLIN,579) 'CDELT2 ', CDelt2*nb 579 FORMAT(A,'= ',F20.8,' / Solar Scale Axis 2.') CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 22 IMBA = IMBA + 20 WRITE(MLIN,5791) 'RX-RY ', deltarad/nb 5791 FORMAT(A,'= ',F20.8,' / Radius x minus Radius y.') CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 23 IMBA = IMBA + 20 WRITE(MLIN,5792) 'GRADXP ', xpslopeave 5792 FORMAT(A,'= ',F20.8,' / Average Slope X plus. ') CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 24 IMBA = IMBA + 20 WRITE(MLIN,5793) 'GRADXM ', xmslopeave 5793 FORMAT(A,'= ',F20.8,' / Average Slope X minus. ') CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 25 IMBA = IMBA + 20 WRITE(MLIN,5794) 'GRADYP ', ypslopeave 5794 FORMAT(A,'= ',F20.8,' / Average Slope Y plus. ') CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 26 IMBA = IMBA + 20 WRITE(MLIN,5795) 'GRADYM ', ymslopeave 5795 FORMAT(A,'= ',F20.8,' / Average Slope Y minus. ') CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 27 IMBA = IMBA + 20 WRITE(MLIN,5796) 'SCSLOPXP', scatxpslope0 5796 FORMAT(A,'= ',F20.8,' / Scattered Light Slope X plus.') CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 28 IMBA = IMBA + 20 WRITE(MLIN,5797) 'SCSLOPXM', scatxmslope0 5797 FORMAT(A,'= ',F20.8,' / Scattered Light Slope X minus.') CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 29 IMBA = IMBA + 20 WRITE(MLIN,5798) 'SCSLOPYP', scatypslope0 5798 FORMAT(A,'= ',F20.8,' / Scattered Light Slope Y plus.') CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 30 IMBA = IMBA + 20 WRITE(MLIN,5799) 'SCSLOPYM', scatymslope0 5799 FORMAT(A,'= ',F20.8,' / Scattered Light Slope Y minus.') CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 31 IMBA = IMBA + 20 do ikw=1, 3 write(mlin,581) mtifkwd(ikw),'= ',mtifkwd_value(ikw) 581 format(A8,A2,A70) call tbyte(bufo(imba),mlin,nca) imba = imba +20 enddo IREC = 1 select case(ioutloc) case(1) WRITE(3,REC=IREC) BUFO case(2) write(13,rec=irec) bufo end select IREC = IREC + 1 DO I=1,720 BUFO(I) = ' ' enddo IMBA = 1 do ikw=4, 38 write(mlin,581) mtifkwd(ikw),'= ',mtifkwd_value(ikw) call tbyte(bufo(imba),mlin,nca) imba = imba +20 enddo WRITE(MLIN,596) 'NPIXDUST',ireplace, 1 'Number of dust and pit pixels replaced.' 596 FORMAT(A,'= ',I20,' / ',A) CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 33 IMBA = IMBA + 20 select case(ioutloc) case(1) WRITE(3,REC=IREC) BUFO case(2) write(13,rec=irec) bufo end select IREC = IREC + 1 DO 595 I=1,720 BUFO(I) = ' ' 595 CONTINUE IMBA = 1 597 FORMAT(A,'= ',f20.7,' / ',A) 598 FORMAT(A,'= ',f20.2,' / ',A) WRITE(MLIN,598) 'RADIUS_X',1.0/cdelt1/nb,'X-direction Radius' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 35 IMBA = IMBA + 20 WRITE(MLIN,598) 'RADIUS_Y',1.0/cdelt2/nb,'Y-direction Radius' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 35 IMBA = IMBA + 20 WRITE(MLIN,597) 'AVE_INTS',valave,'Average Intensity Value' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 35 IMBA = IMBA + 20 WRITE(MLIN,597) 'AVE_IRMS',amsave, 1 'Ave. RMS Intensity Variation over 5x5 Arrays' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 35 IMBA = IMBA + 20 WRITE(MLIN,597) 'SUNEAR_F',sun_earth_fac, 1 'Sun-Earth Distance Factor' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 34 IMBA = IMBA + 20 write(mlin,596) 'RADSCORE',idemerit, 1 'Radius Solution Demerits, < 6 is good' call tbyte(bufo(imba), mlin, nca) imba = imba + 20 WRITE(MLIN,30) 'END ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 36 IMBA = IMBA + 20 select case (ioutloc) case(1) WRITE(3,REC=IREC) BUFO case(2) write(13,rec=irec) bufo end select IREC = IREC + 1 ndxob = 0 DO 670 J=1,NAX2A/nb ILI = (J - 1)*NAX1A/nb do 660 i = 1,nax1a/nb if (nb.gt.1) then fpic(ili+i) = 0.0 do 656 jb = 1,nb do 654 ib = 1,nb fpic(ili+i) = fpic(ili+i) + arr(nb*(i-1)+ib,nb*(j-1)+jb) 654 continue 656 continue ipic(ili+i)= nint(fpic(ili+i)/nb**2) else IPIC(ILI + I)=arr(i,j) endif ndxob = ndxob + 1 660 CONTINUE 670 CONTINUE C DUMP IPIC ARRAY ON DISK IN FITS FORMAT. 700 K = 1 DO 720 L=1,NDXOb IF(K.LE.1440) GO TO 710 CALL RESWAB(IBFO, 1440) select case(ioutloc) case(1) WRITE(3,REC=IREC) IBFO case(2) write(13,rec=irec) ibfo end select IREC = IREC + 1 K = 1 710 IBFO(K) = IPIC(L) K = K + 1 720 CONTINUE IF(K.GT.1440) GO TO 740 DO 730 I=K,1440 IBFO(I) = NUL 730 CONTINUE CALL RESWAB(IBFO, 1440) select case(ioutloc) case(1) WRITE(3,REC=IREC) IBFO case(2) write(13,rec=irec) ibfo end select 740 continue 745 continue CLOSE(UNIT=3) close(unit=13) C WRITE dusty FITS HEADER. c nb=5 DO I=1,720 BUFO(I) = ' ' end do IMBA = 1 WRITE(MLIN,530) 'SIMPLE ', 'T', 'FITS STANDARD' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 1 IMBA = IMBA + 20 WRITE(MLIN,540) 'BITPIX ', IBITPIX, 'FITS BITS/PIXEL' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 2 IMBA = IMBA + 20 WRITE(MLIN,540) 'NAXIS ', 2, 'NUMBER OF AXES' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 3 IMBA = IMBA + 20 WRITE(MLIN,540) 'NAXIS1 ', NAX1A/nb, 'X AXIS NO. PIX' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 4 IMBA = IMBA + 20 WRITE(MLIN,540) 'NAXIS2 ', NAX2A/nb, 'Y AXIS NO. PIX' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 5 IMBA = IMBA + 20 WRITE(MLIN,540) 'BLANK ', -32768, 'ARRAY NULL VALUE' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 6 IMBA = IMBA + 20 write(mlin,545) iyearrun,imonthrun,idayrun, 1 ihourrun,iminrun,isecrun CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 7 IMBA = IMBA + 20 write(mlin,746) 746 format( 1 'HISTORY Raw image with dust and emulsion pits remaining on the 2scan. ') CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 7 IMBA = IMBA + 20 WRITE(MLIN,550) 'BSCALE ', 1.0, 'REAL = DATA*BSCALE + BZERO' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 7 IMBA = IMBA + 20 WRITE(MLIN,550) 'BZERO ', 0.0, ' ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 8 IMBA = IMBA + 20 WRITE(MLIN,530) 'BUNIT ', '''INTENS. WITH DUST'' ', 1 'ARBITRARY UNITS' imba_quant = imba CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 9 IMBA = IMBA + 20 WRITE(MLIN,530) 'ORIGIN ', '''MT. WILSON'' ', ' ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 10 IMBA = IMBA + 20 WRITE(MLIN,530) 'TELESCOP', '''60 FT'' ', ' ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 11 IMBA = IMBA + 20 imba_obj = imba WRITE(MLIN,530) 'OBJECT ', '''SUN'' ', ' ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 12 IMBA = IMBA + 20 IF(WAVEL.EQ.0.0) THEN WRITE(MLIN,530) 'DATA-TYP', '''DIRECT'' ', ' ' ELSE WRITE(MLIN,530) 'DATA-TYP', '''CA SPECTROHELIOGRM''', ' ' ENDIF CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 13 IMBA = IMBA + 20 WRITE(MLIN,560) 'DATE-OBS', IYR, IMO, IDA, 'YEAR-MO-DA' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 14 IMBA = IMBA + 20 WRITE(MLIN,550) 'TAVG-UT ', TAV, 'AVERAGE TIME, U.T.' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 15 IMBA = IMBA + 20 IF(WAVEL.NE.0.0) THEN WRITE(MLIN,550) 'LAMBDA ', WAVEL, 'WAVELENGTH (ANGSTROMS)' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 16 IMBA = IMBA + 20 ENDIF WRITE(MLIN,570) 'SCALPIX1', SclPIX1/nb, 'X', MRESU CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 17 IMBA = IMBA + 20 WRITE(MLIN,570) 'SCALPIX2', SclPIX2/nb, 'Y', MRESU CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 18 IMBA = IMBA + 20 WRITE(MLIN,575) 'CRPIX1 ', CRPIX1/nb CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 19 IMBA = IMBA + 20 WRITE(MLIN,576) 'CRPIX2 ', CRPIX2/nb CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 20 IMBA = IMBA + 20 WRITE(MLIN,577) 'CDELT1 ', CDelt1*nb CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 21 IMBA = IMBA + 20 WRITE(MLIN,579) 'CDELT2 ', CDelt2*nb CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 22 IMBA = IMBA + 20 WRITE(MLIN,5791) 'RX-RY ', deltarad CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 23 IMBA = IMBA + 20 WRITE(MLIN,5792) 'GRADXP ', xpslopeave CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 24 IMBA = IMBA + 20 WRITE(MLIN,5793) 'GRADXM ', xmslopeave CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 25 IMBA = IMBA + 20 WRITE(MLIN,5794) 'GRADYP ', ypslopeave CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 26 IMBA = IMBA + 20 WRITE(MLIN,5795) 'GRADYM ', ymslopeave CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 27 IMBA = IMBA + 20 WRITE(MLIN,5796) 'SCSLOPXP', scatxpslope0 CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 28 IMBA = IMBA + 20 WRITE(MLIN,5797) 'SCSLOPXM', scatxmslope0 CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 29 IMBA = IMBA + 20 WRITE(MLIN,5798) 'SCSLOPYP', scatypslope0 CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 30 IMBA = IMBA + 20 WRITE(MLIN,5799) 'SCSLOPYM', scatymslope0 CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 31 IMBA = IMBA + 20 do ikw=1, 3 write(mlin,581) mtifkwd(ikw),'= ',mtifkwd_value(ikw) call tbyte(bufo(imba),mlin,nca) imba = imba +20 enddo IREC = 1 WRITE(4,REC=IREC) BUFO IREC = IREC + 1 DO I=1,720 BUFO(I) = ' ' enddo IMBA = 1 do ikw=4, 38 write(mlin,581) mtifkwd(ikw),'= ',mtifkwd_value(ikw) call tbyte(bufo(imba),mlin,nca) imba = imba +20 enddo WRITE(MLIN,596) 'NPIXDUST',ireplace CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 33 IMBA = IMBA + 20 WRITE(4,REC=IREC) BUFO IREC = IREC + 1 DO 753 I=1,720 BUFO(I) = ' ' 753 CONTINUE IMBA = 1 WRITE(MLIN,598) 'RADIUS_X',1.0/cdelt1 CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 35 IMBA = IMBA + 20 WRITE(MLIN,598) 'RADIUS_Y',1.0/cdelt2 CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 35 IMBA = IMBA + 20 WRITE(MLIN,597) 'AVE_INTS',valave CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 35 IMBA = IMBA + 20 WRITE(MLIN,597) 'AVE_IRMS',amsave CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 35 IMBA = IMBA + 20 WRITE(MLIN,597) 'SUNEAR_F',sun_earth_fac CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 34 IMBA = IMBA + 20 write(mlin,596) 'RADSCORE',idemerit call tbyte(bufo(imba), mlin, nca) imba = imba + 20 WRITE(MLIN,30) 'END ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 23 IMBA = IMBA + 20 c NLINE = (IMBA - 1)/20 c J = 1 c DO I=1,NLINE c TYPE 580, (BUFO(K),K=J,J+19) c J = J + 20 c end do WRITE(4,REC=IREC) BUFO IREC = IREC + 1 ndxob = 0 DO 770 J=1,NAX2A/nb ILI = (J - 1)*(NAX1A/nb) DO 760 I=1,NAX1A/nb fpic(ili+i) = 0.0 do 756 jb=1,nb do 754 ib=1,nb fpic(ili+i) = fpic(ili+i) + arr0(nb*(i-1)+ib,nb*(j-1)+jb) 754 continue 756 continue IPIC(ILI + I)=nint(fpic(ili+i)/nb**2) ndxob = ndxob + 1 c IPIC(ILI + I)=arr0(i,j) 760 CONTINUE 770 CONTINUE C DUMP IPIC ARRAY ON DISK IN FITS FORMAT. 800 K = 1 DO 820 L=1,NDXOb IF(K.LE.1440) GO TO 810 CALL RESWAB(IBFO, 1440) WRITE(4,REC=IREC) IBFO IREC = IREC + 1 K = 1 810 IBFO(K) = IPIC(L) K = K + 1 820 CONTINUE IF(K.GT.1440) GO TO 840 DO 830 I=K,1440 IBFO(I) = NUL 830 CONTINUE CALL RESWAB(IBFO, 1440) WRITE(4,REC=IREC) IBFO 840 CLOSE(UNIT=4) c c C WRITE rms FITS HEADER. c nb=5 DO I=1,720 BUFO(I) = ' ' end do IMBA = 1 WRITE(MLIN,530) 'SIMPLE ', 'T', 'FITS STANDARD' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 1 IMBA = IMBA + 20 WRITE(MLIN,540) 'BITPIX ', IBITPIX, 'FITS BITS/PIXEL' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 2 IMBA = IMBA + 20 WRITE(MLIN,540) 'NAXIS ', 2, 'NUMBER OF AXES' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 3 IMBA = IMBA + 20 WRITE(MLIN,540) 'NAXIS1 ', NAX1A/nb, 'X AXIS NO. PIX' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 4 IMBA = IMBA + 20 WRITE(MLIN,540) 'NAXIS2 ', NAX2A/nb, 'Y AXIS NO. PIX' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 5 IMBA = IMBA + 20 WRITE(MLIN,540) 'BLANK ', -32768, 'ARRAY NULL VALUE' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 6 IMBA = IMBA + 20 WRITE(MLIN,550) 'BSCALE ', 1.0, 'REAL = DATA*BSCALE + BZERO' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 7 IMBA = IMBA + 20 WRITE(MLIN,550) 'BZERO ', 0.0, ' ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 8 IMBA = IMBA + 20 WRITE(MLIN,530) 'BUNIT ', '''RMS INTENSITY DEV.''', 1 'ARBITRARY UNITS' imba_quant = imba CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 9 IMBA = IMBA + 20 WRITE(MLIN,530) 'ORIGIN ', '''MT. WILSON'' ', ' ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 10 IMBA = IMBA + 20 WRITE(MLIN,530) 'TELESCOP', '''60 FT'' ', ' ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 11 IMBA = IMBA + 20 imba_obj = imba WRITE(MLIN,530) 'OBJECT ', '''SUN'' ', ' ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 12 IMBA = IMBA + 20 IF(WAVEL.EQ.0.0) THEN WRITE(MLIN,530) 'DATA-TYP', '''DIRECT'' ', ' ' ELSE WRITE(MLIN,530) 'DATA-TYP', '''CA SPECTROHELIOGRM''', ' ' ENDIF CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 13 IMBA = IMBA + 20 WRITE(MLIN,560) 'DATE-OBS', IYR, IMO, IDA, 'YEAR-MO-DA' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 14 IMBA = IMBA + 20 WRITE(MLIN,550) 'TAVG-UT ', TAV, 'AVERAGE TIME, U.T.' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 15 IMBA = IMBA + 20 IF(WAVEL.NE.0.0) THEN WRITE(MLIN,550) 'LAMBDA ', WAVEL, 'WAVELENGTH (ANGSTROMS)' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 16 IMBA = IMBA + 20 ENDIF WRITE(MLIN,570) 'SCALPIX1', SclPIX1/nb, 'X', MRESU CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 17 IMBA = IMBA + 20 WRITE(MLIN,570) 'SCALPIX2', SclPIX2/nb, 'Y', MRESU CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 18 IMBA = IMBA + 20 WRITE(MLIN,575) 'CRPIX1 ', CRPIX1/nb CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 19 IMBA = IMBA + 20 WRITE(MLIN,576) 'CRPIX2 ', CRPIX2/nb CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 20 IMBA = IMBA + 20 WRITE(MLIN,577) 'CDELT1 ', CDelt1*nb CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 21 IMBA = IMBA + 20 WRITE(MLIN,579) 'CDELT2 ', CDelt2*nb CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 22 IMBA = IMBA + 20 WRITE(MLIN,30) 'END ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 23 IMBA = IMBA + 20 c NLINE = (IMBA - 1)/20 c J = 1 c DO I=1,NLINE c TYPE 580, (BUFO(K),K=J,J+19) c J = J + 20 c end do IREC = 1 WRITE(2,REC=IREC) BUFO IREC = IREC + 1 ndxob = 0 DO 1070 J=1,NAX2A/nb ILI = (J - 1)*(NAX1A/nb) DO 1060 I=1,NAX1A/nb fpic(ili+i) = 0.0 do 1056 jb=1,nb do 1054 ib=1,nb fpic(ili+i) = fpic(ili+i) + avearr(nb*(i-1)+ib,nb*(j-1)+jb) 1054 continue 1056 continue IPIC(ILI + I)=nint(fpic(ili+i)/nb**2) ndxob = ndxob + 1 c IPIC(ILI + I)=avearr(i,j) 1060 CONTINUE 1070 CONTINUE C DUMP IPIC ARRAY ON DISK IN FITS FORMAT. K = 1 DO 1120 L=1,NDXOb IF(K.LE.1440) GO TO 1110 CALL RESWAB(IBFO, 1440) WRITE(2,REC=IREC) IBFO IREC = IREC + 1 K = 1 1110 IBFO(K) = IPIC(L) K = K + 1 1120 CONTINUE IF(K.GT.1440) GO TO 1140 DO 1130 I=K,1440 IBFO(I) = NUL 1130 CONTINUE CALL RESWAB(IBFO, 1440) WRITE(2,REC=IREC) IBFO 1140 CLOSE(UNIT=2) c c C WRITE slope FITS HEADER. c DO I=1,720 BUFO(I) = ' ' end do IMBA = 1 WRITE(MLIN,530) 'SIMPLE ', 'T', 'FITS STANDARD' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 1 IMBA = IMBA + 20 WRITE(MLIN,540) 'BITPIX ', IBITPIX, 'FITS BITS/PIXEL' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 2 IMBA = IMBA + 20 WRITE(MLIN,540) 'NAXIS ', 2, 'NUMBER OF AXES' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 3 IMBA = IMBA + 20 WRITE(MLIN,540) 'NAXIS1 ', NAX1A/nb, 'X AXIS NO. PIX' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 4 IMBA = IMBA + 20 WRITE(MLIN,540) 'NAXIS2 ', NAX2A/nb, 'Y AXIS NO. PIX' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 5 IMBA = IMBA + 20 WRITE(MLIN,540) 'BLANK ', -32768, 'ARRAY NULL VALUE' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 6 IMBA = IMBA + 20 WRITE(MLIN,550) 'BSCALE ', 1.0, 'REAL = DATA*BSCALE + BZERO' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 7 IMBA = IMBA + 20 WRITE(MLIN,550) 'BZERO ', 0.0, ' ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 8 IMBA = IMBA + 20 WRITE(MLIN,530) 'BUNIT ', '''INTENSITY GRAD .''', 1 'ARBITRARY UNITS' imba_quant = imba CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 9 IMBA = IMBA + 20 WRITE(MLIN,530) 'ORIGIN ', '''MT. WILSON'' ', ' ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 10 IMBA = IMBA + 20 WRITE(MLIN,530) 'TELESCOP', '''60 FT'' ', ' ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 11 IMBA = IMBA + 20 imba_obj = imba WRITE(MLIN,530) 'OBJECT ', '''SUN'' ', ' ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 12 IMBA = IMBA + 20 IF(WAVEL.EQ.0.0) THEN WRITE(MLIN,530) 'DATA-TYP', '''DIRECT'' ', ' ' ELSE WRITE(MLIN,530) 'DATA-TYP', '''CA SPECTROHELIOGRM''', ' ' ENDIF CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 13 IMBA = IMBA + 20 WRITE(MLIN,560) 'DATE-OBS', IYR, IMO, IDA, 'YEAR-MO-DA' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 14 IMBA = IMBA + 20 WRITE(MLIN,550) 'TAVG-UT ', TAV, 'AVERAGE TIME, U.T.' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 15 IMBA = IMBA + 20 IF(WAVEL.NE.0.0) THEN WRITE(MLIN,550) 'LAMBDA ', WAVEL, 'WAVELENGTH (ANGSTROMS)' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 16 IMBA = IMBA + 20 ENDIF WRITE(MLIN,570) 'SCALPIX1', SclPIX1/nb, 'X', MRESU CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 17 IMBA = IMBA + 20 WRITE(MLIN,570) 'SCALPIX2', SclPIX2/nb, 'Y', MRESU CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 18 IMBA = IMBA + 20 WRITE(MLIN,575) 'CRPIX1 ', CRPIX1/nb CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 19 IMBA = IMBA + 20 WRITE(MLIN,576) 'CRPIX2 ', CRPIX2/nb CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 20 IMBA = IMBA + 20 WRITE(MLIN,577) 'CDELT1 ', CDelt1*nb CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 21 IMBA = IMBA + 20 WRITE(MLIN,579) 'CDELT2 ', CDelt2*nb CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 22 IMBA = IMBA + 20 WRITE(MLIN,30) 'END ' CALL TBYTE(BUFO(IMBA), MLIN, NCA) ! card 23 IMBA = IMBA + 20 c NLINE = (IMBA - 1)/20 c J = 1 c DO I=1,NLINE c TYPE 580, (BUFO(K),K=J,J+19) c J = J + 20 c end do IREC = 1 WRITE(7,REC=IREC) BUFO IREC = IREC + 1 ndxob = 0 DO 1170 J=1,NAX2A/nb ILI = (J - 1)*(NAX1A/nb) DO 1160 I=1,NAX1A/nb c IPIC(ILI + I)=difarr(i,j) fpic(ili+i) = 0.0 do 1156 jb=1,nb do 1154 ib=1,nb fpic(ili+i) = fpic(ili+i) + difarr(nb*(i-1)+ib,nb*(j-1)+jb) 1154 continue 1156 continue IPIC(ILI + I)=nint(fpic(ili+i)/nb**2) ndxob = ndxob + 1 c IPIC(ILI + I)=avearr(i,j) 1160 CONTINUE 1170 CONTINUE C DUMP IPIC ARRAY ON DISK IN FITS FORMAT. K = 1 DO 1220 L=1,NDXOb IF(K.LE.1440) GO TO 1210 CALL RESWAB(IBFO, 1440) WRITE(7,REC=IREC) IBFO IREC = IREC + 1 K = 1 1210 IBFO(K) = IPIC(L) K = K + 1 1220 CONTINUE IF(K.GT.1440) GO TO 1240 DO 1230 I=K,1440 IBFO(I) = NUL 1230 CONTINUE CALL RESWAB(IBFO, 1440) WRITE(7,REC=IREC) IBFO 1240 CLOSE(UNIT=7) c c TYPE *, 'NO. REC', IREC c c ! Free the memory allocation goto 35 880 continue close(unit=11) c close(unit=12) if (allocated) then deallocate(ARR) deallocate(avearr) deallocate(difarr) deallocate(mask) deallocate(arr0) deallocate(xave) deallocate(yave) deallocate(xslopeave) deallocate(yslopeave) deallocate(slopeave) endif GO TO 950 890 TYPE *, 'BAD TIF FILE' 950 STOP END SUBROUTINE GETIFEL(IBADR, IARR, NBYT, ISWBF, LUN) C C GETS A TIFF FILE ELEMENT. C IBADR - BYTE ADDRESS OF ELEMENT REFERENCED TO 1ST BYTE IN FILE, 0 . C IARR - INTEGER*2 ARRAY TO PUT ELEMENT IN. C NBYT - NO. OF BYTES TO TRANSFER. C ISWBF - SWAP BYTES FLAG. =0 DON'T. =1, DO. C LUN - LOGICAL UNIT TO HAVE A DIRECT ACCESS UNFORMATTED TIF C FILE OPENED WITH 512 BYTE RECL (RECL=128 ON DEC). C INTEGER*2 IAR(256), IARR((NBYT+1)/2) DATA ICREC/ 0/ IREC = IBADR/512 + 1 IF(IREC.EQ.ICREC) GO TO 20 C NEED TO READ. READ(LUN,REC=IREC,ERR=50) IAR IF(ISWBF.NE.0) CALL SWAB(IAR, 512) ICREC = IREC 20 IWD = MOD(IBADR,512)/2 + 1 NWT = (NBYT + 1)/2 DO 40 I=1,NWT IF(IWD.LE.256) GO TO 30 C SPANS END OF RECORD, NEED TO READ AGAIN. IREC = IREC + 1 READ(LUN,REC=IREC,ERR=50) IAR IF(ISWBF.NE.0) CALL SWAB(IAR, 512) ICREC = IREC IWD = 1 30 IARR(I) = IAR(IWD) IWD = IWD + 1 40 CONTINUE 50 CONTINUE RETURN END SUBROUTINE IFDNTRY(IFDA, ISWBF, LUN, ITAG, ITYP, ICOUNT, IDATA, 1 RDATA, IVAF) C C RETURNS AN IMAGE FILE DIRECTORY ENTRY (IFD) FROM A TIF FILE. C C IFDA - IFD ADDRESS, MUST BE A VARIABLE, UPDATED TO NEXT IFD ADDRESS. C ISWBF - SWAP BYTES FLAG. =0 DON'T. =1, DO. C LUN - LOGICAL UNIT TO HAVE A DIRECT ACCESS UNFORMATTED TIF C FILE OPENED WITH 512 BYTE RECL (RECL=128 ON DEC). C ITAG - THE TAG THAT IDENTIFIES THE FIELD. C ITYP - FIELD TYPE: =1, BYTE. =2, ASCII. =3, INTEGER*2. C =4, INTEGER*4. =5, RATIONAL (SEE RDATA AND IVAF, BELOW). C ICOUNT - NUMBER OF VALUES. C IDATA - EITHER THE VALUE OR THE ADDRESS OF THE VALUE. C RDATA - IF ITYP=5 AND ICOUNT=1, A FLOATING POINT VALUE. C IVAF - FLAG. =0, A VALUE IS RETURNED IN IDATA OR RDATA. C =1, THE VALUE(S) ARE TOO MANY OR TOO LONG TO FIT IN C A 4 BYTE WORD (IDATA OR RDATA), IDATA POINTS TO THE C ARRAY. IN THE CASE OF ITYP=5, EACH VALUE IS A PAIR OF C INTEGER*4 NUMBERS: 1ST NUMBER / 2ND NUMBER = FLOATING C POINT VALUE. C INTEGER*2 IARR(6) RDATA = 0. CALL GETIFEL(IFDA, IARR, 12, ISWBF, LUN) ITAG = IFOUR(IARR(1)) ITYP = IARR(2) ICOUNT = I2TOI4(IARR(3)) IVAF = 1 IF((ITYP.LT.3.AND.ICOUNT.LE.4).OR. 1 (ITYP.EQ.3.AND.ICOUNT.LE.2).OR. 2 (ITYP.EQ.4.AND.ICOUNT.LE.1)) IVAF = 0 IDATA = I2TOI4(IARR(5)) IF(ICOUNT.GT.1.OR.ITYP.NE.5) GO TO 30 CALL GETIFEL(IDATA, IARR, 8, ISWBF, LUN) RDATA = FLOAT(I2TOI4(IARR(1)))/FLOAT(I2TOI4(IARR(3))) IVAF = 0 30 IFDA = IFDA + 12 RETURN END FUNCTION I2TOI4(I2A) C MAKE 2 I*2 WORDS INTO VAX I*4 WORD. INTEGER*2 I2A(2), IREF(2) EQUIVALENCE (IREF,I4W) IREF(1) = I2A(1) IREF(2) = I2A(2) I2TOI4 = I4W RETURN END SUBROUTINE RTJFIL(MARR, MFIL, NCFW) C C RIGHT JUSTIFIES THE CHARACTERS IN MARR(1:NCFW) AND FILLS WITH C THE CHARACTER IN MFIL. C CHARACTER*(*) MARR, MFIL CHARACTER*1 MF2 MF2 = MFIL(1:1) NCIA = LEN_TRIM(MARR(1:NCFW)) NCT = NCFW + 1 - NCIA IF(NCT.LT.1) GO TO 20 MARR(NCT:NCFW) = MARR(1:NCIA) DO 10 I=1,NCT-1 MARR(I:I) = MF2 10 CONTINUE 20 RETURN END SUBROUTINE WREBIN(AL, NXA, NYA, INIT, SL, WT, NXS, NYS, NSO) C C REBINS A 2 DIMENSIONAL ARRAY TO FORM ANOTHER 2 DIMENSIONAL ARRAY C OF A DIFFERENT SIZE. THIS IS DONE BY WEIGHTED SUMS. THE INPUT C ARRAY AL, IS INPUT 1 LINE AT A TIME, THE OUTPUT ARRAY, SL IS C OUTPUT 1 LINE AT A TIME. IF OUTPUT INTEGER NSO IS 0, WREBIN C NEEDS ANOTHER LINE INPUT (REAL VARIABLE AL) BUT THERE IS NO LINE C OUTPUT (REAL VARIABLE SL). IF NSO IS 1 OR GREATER, WREBIN HAS C AN OUTPUT LINE IN SL THAT NEEDS TO BE USED NSO TIMES BUT REQUIRES C THE SAME LINE IN AL AS BEFORE. C AL - INPUT LINE OF A 2 DIMENSIONAL ARRAY C NXA - NO. OF VALUES IN AL. C NYA - NO OF LINES IN THE INPUT ARRAY C INIT - MUST BE A VARIABLE AND MUST HAVE A NON-ZERO VALUE FOR C THE FIRST CALL TO WREBIN. C SL - RE-BINNED OUTPUT LINE. C WT - NEEDED BY WREBIN TO STORE WEIGHTS, SAME SIZE AS SL C NXS - NO. OF VALUES TO STORE IN ARRAYS SL AND WT. C NYS - NO OF LINES IN RE-BINNED ARRAY. C NSO - OUTPUT FLAG: =0, NEW INPUT LINE NEEDED, NO OUTPUT LINE. C >0, OUTPUT LINE NEEDS TO BE USED NSO TIMES, LEAVE OLD C INPUT LINE. DON'T CHANGE THE VALUE OF NSO. C DIMENSION AL(NXA+1), SL(NXS+1), WT(NXS+1) IF(INIT.EQ.0) GO TO 20 C INITIALIZE. ANUL = RNUL(0) INIT = 0 XAI = FLOAT(NXS)/FLOAT(NXA) XSI = 1. YSMAX = NYS YAI = YSMAX/FLOAT(NYA) YSI = 1. YAO = 0. YA = YAI YSO = 0. YS = YSI YWT = YSI IF(YAI.LT.YSI) YWT = YAI XWT = XSI IF(XAI.LT.XSI) XWT = XAI HAPX = XWT*YWT/2. NSO = 1 20 IF(NSO.EQ.0) GO TO 24 DO 22 I=1,NXS SL(I) = 0. WT(I) = 0. 22 CONTINUE 24 IF(YWT.EQ.0..OR.YS.GT.YSMAX) GO TO 160 I = 1 J = 1 XAO = 0. XA = XAI XSO = 0. XS = XSI IF(XAI.LT.XSI) GO TO 100 C SCREEN (OUTPUT) PIXELS SMALLER IN X. 30 DIF = XA - XSO IF(DIF.GE.XSI) GO TO 40 ANT = AL(J) IF(ANT.EQ.ANUL) GO TO 32 WTT = DIF*YWT WT(I) = WT(I) + WTT SL(I) = SL(I) + ANT*WTT 32 J = J + 1 XAO = XA XA = XA + XAI DIF = XS - XAO WTT = DIF*YWT GO TO 50 40 WTT = XSI*YWT 50 ANT = AL(J) IF(ANT.EQ.ANUL) GO TO 52 WT(I) = WT(I) + WTT SL(I) = SL(I) + ANT*WTT 52 I = I + 1 IF(I.GT.NXS) GO TO 130 XSO = XS XS = XS + XSI GO TO 30 C SCREEN PIXELS LARGER IN X. 100 DIF = XS - XAO IF(DIF.GE.XAI) GO TO 110 ANT = AL(J) IF(ANT.EQ.ANUL) GO TO 102 WTT = DIF*YWT WT(I) = WT(I) + WTT SL(I) = SL(I) + ANT*WTT 102 I = I + 1 IF(I.GT.NXS) GO TO 130 XSO = XS XS = XS + XSI DIF = XA - XSO WTT = DIF*YWT GO TO 120 110 WTT = XAI*YWT 120 ANT = AL(J) IF(ANT.EQ.ANUL) GO TO 122 WT(I) = WT(I) + WTT SL(I) = SL(I) + ANT*WTT 122 J = J + 1 XAO = XA XA = XA + XAI GO TO 100 C DONE WITH A SCAN. 130 IF(YS.GT.YA.AND.YS.LT.YSMAX) GO TO 160 DO 140 I=1,NXS ANT = WT(I) IF(ANT.GE.HAPX) GO TO 132 SL(I) = ANUL GO TO 140 132 SL(I) = SL(I)/ANT 140 CONTINUE RNSO = (YA - YSO)/YSI IF(RNSO.GT.YAI.OR.RNSO.LT.1.) RNSO = 1. NSO = RNSO YS = YS + YSI*FLOAT(NSO) YSO = YS - YSI YWT = YA - YSO IF(YWT.GT.YSI) YWT = YSI GO TO 180 160 NSO = 0 YAO = YA YA = YA + YAI YWT = YS - YAO IF(YWT.GT.YAI) YWT = YAI 180 RETURN END SUBROUTINE HISTIN(AR, NBIN, VMIN, VMAX, ISUB) C C INITIALIZES THE HISTOGRAM ROUTINE HISTDO. C AR - ARRAY TO BIN IN. C NBIN - NO. VALUES IN AR. C VMIN, VMAX - MINIMUM AND MAXIMUM VALUES TO CONSIDER. C ISUB - CALL SUBSCRIPT, NORMALLY 1 UNLESS MULTIPLE C HISTOGRAMS ARE COMPUTED IN THE SAME LOOP. C MAX 50 . C DIMENSION EM(50), BE(50), NBINS(50), AR(NBIN), VAL(NV) SAVE EM, BE, NBINS, ANUL ANUL = RNUL(0) NBINS(ISUB) = NBIN DO 10 I=1,NBIN AR(I) = 0. 10 CONTINUE EM(ISUB) = FLOAT(NBIN)/(VMAX - VMIN) BE(ISUB) = -VMIN*EM(ISUB) GO TO 30 ENTRY HISTDO(VAL, NV, AR, NBIN, ISUB) C C DOES A HISTOGRAM FOR ARRAY VAL. MULTIPLE CALLS WILL ADD TO C THE HISTOGRAM. C VAL - ARRAY OF VALUES TO USE FOR HISTOGRAM. C NV - NO OF VAL'S FOR THIS CALL. C AR - HISTOGRAM ARRAY. C ISUB - SEE ISUB ABOVE. C EML = EM(ISUB) BEL = BE(ISUB) NBINL = NBINS(ISUB) DO 20 I=1,NV ANT = VAL(I) IF(ANT.EQ.ANUL) GO TO 20 J = EML*ANT + BEL + 1.0 IF(J.LT.1) J = 1 IF(J.GT.NBINL) J = NBINL AR(J) = AR(J) + 1. 20 CONTINUE 30 RETURN END INTEGER*4 FUNCTION IFOUR(IDUM2) C RETURNS I*4 FROM UNSIGNED I*2 INTEGER*2 ITEM(2), IDUM2 INTEGER*4 JTEM EQUIVALENCE (ITEM,JTEM) ITEM(1) = IDUM2 IFOUR = JTEM RETURN END FUNCTION LEN_TRIM(MCHR) C C RETURNS THE LENGTH OF A CHARACTER ARRAY MINUS THE TRAILING BLANKS. C CHARACTER*(*) MCHR NCA = LEN(MCHR) DO 10 I=NCA,1,-1 IF(MCHR(I:I).EQ.' ') GO TO 10 LEN_TRIM = I GO TO 20 10 CONTINUE LEN_TRIM = 0 20 RETURN END SUBROUTINE RESWAB(LAR, NI2) C C SWAPS EACH PAIR OF BYTES IN 4 BYTE WORDS SO THAT THE 4,3,2,1 (DEC) ORDER C BECOMES 3,4,1,2 (SUN) OR VICE-VERSA. NI2 IS THE NO. OF I*2 WORDS IN THE C I*4 ARRAY C BYTE MAR(4), MR1, MR2, MR3, MR4 DIMENSION LAR((NI2+1)/2) EQUIVALENCE (MAR,LTM), (MR1,MAR(1)), (MR2,MAR(2)), (MR3,MAR(3)), 1 (MR4,MAR(4)) NI4 = NI2/2 DO 10 I=1,NI4 LTM = LAR(I) INT = MR1 !TRADE MR1 AND MR2. MR1 = MR2 MR2 = INT INT = MR3 !TRADE MR3 AND MR4. MR3 = MR4 MR4 = INT LAR(I) = LTM 10 CONTINUE IF(MOD(NI2,2).NE.0) THEN LTM = LAR(NI4+1) INT = MR1 MR1 = MR2 MR2 = INT LAR(NI4+1) = LTM ENDIF RETURN END SUBROUTINE RENOSB(LUN, RNOA, NNR, NNMAX, MDELA) C C RETURNS REALS FROM NUMBERS READ FROM LUN. C RNOA IS THE REAL ARRAY TO RECEIVE DATA. C NNR IS THE NO. OF REALS RETURNED. C NNMAX IS THE MAXIMUM NO. OF REALS TO RETURN. C MDELA IS A CHARACTER ARRAY OF DELIMITERS. C A MAXIMUM OF 132 CHARACTERS, READ FROM LUN, ARE PROCESSED PER CALL. C IMBEDDED + OR - IS TREATED FIRST AS A DELIMITER AND SECOND C AS THE SIGN OF THE FOLLOWING NO. TWO ADJACENT DELIMITERS OR C A LEADING DELIMITER RETURNS A FLOATING POINT NULL C (SEE FUNCTION RNUL). C CHARACTER*132 MWRK CHARACTER*(*) MDELA CHARACTER*20 MTMP CHARACTER*1 MCH DIMENSION RNOA(NNMAX) NNR = 0 NDEL = LEN(MDELA) READ(LUN,10,END=110) MWRK 10 FORMAT(A) NC = LEN_TRIM(MWRK) IF(NC.EQ.0) GO TO 110 J = 0 DO 100 I=1,NC INSIGN = 0 MCH = MWRK(I:I) IF(MCH.NE.'+'.AND.MCH.NE.'-') GO TO 20 INSIGN = 1 IF(J.EQ.0) GO TO 40 K = I - 1 IF(MWRK(K:K).EQ.'E'.OR.MWRK(K:K).EQ.'e') GO TO 40 INSIGN = 2 GO TO 50 20 DO 30 L=1,NDEL IF(MCH.EQ.MDELA(L:L)) GO TO 50 30 CONTINUE 40 J = J + 1 MTMP(J:J) = MCH IF(I.LT.NC.OR.INSIGN.EQ.2) GO TO 100 50 IF(MCH.EQ.' '.AND.J.EQ.0) GO TO 100 IF(NNR.LT.NNMAX) GO TO 60 J = 0 GO TO 100 60 NNR = NNR + 1 IF(J.NE.1.OR.INSIGN.NE.1) GO TO 62 J = 0 RNOA(NNR) = RNUL(0) GO TO 100 62 IF(J.GT.0) GO TO 70 RNOA(NNR) = RNUL(0) GO TO 90 70 READ(MTMP(1:J),80) RNOA(NNR) 80 FORMAT(F20.0) 90 J = 0 IF(INSIGN.NE.0) GO TO 40 100 CONTINUE 110 RETURN END FUNCTION RNUL(IDUM) C INTEGER*4 INUL C EQUIVALENCE(INUL,ANUL) C DATA INUL/ ZFFFEFFFF/ RNUL = -1.7014116E+38 RETURN END SUBROUTINE SWAB(MBF,INB) C C THIS SUBROUTINE SWAPS EACH PAIR OF BYTES IN ARRAY MBF. INB IS THE NO. C OF BYTES IN ARRAY MBF TO BE SWAPPED. INB MUST BE EVEN. C LOGICAL*1 MBF(INB),M DO 10 I=2,INB,2 M=MBF(I) J=I-1 MBF(I)=MBF(J) MBF(J)=M 10 CONTINUE RETURN END SUBROUTINE TBYTE(MBYTA, MCHR, NCA) C C TRANSLATES A CHARACTER VARIABLE TO A BYTE ARRAY. C MBYTA - TRANSLATED BYTE ARRAY (MUST HAVE ENOUGH ELEMENTS TO C RECEIVE ALL OF THE CHARACTERS). C MCHR - INPUT CHARACTER VARIABLE. C NCA - RETURNED NO. OF CHARACTERS IN MBYTA. C CHARACTER*(*) MCHR BYTE MBYTA(*) NCA = LEN_TRIM(MCHR) IF(NCA.LE.0) GO TO 20 DO 10 I=1,NCA MBYTA(I) = ICHAR(MCHR(I:I)) 10 CONTINUE 20 RETURN END SUBROUTINE CHARLN(MBYTA, NCBA, NCLEF, MLIN, NCML) C C BREAKS UP A BYTE ARRAY WITH IMBEDDED LINEFEEDS INTO INDIVIDUAL LINES. C EACH LINE IS RETURNED AS A CHARACTER ARRAY, 1 FOR EACH CALL. C C MBYTA - INPUT BYTE ARRAY. C NCBA - NO. OF CHARACTERS IN MBYTA (INPUT). C NCLEF - NO. OF CHARACTERS LEFT IN MBYTA, MUST BE INITIALIZED TO C NCBA BEFORE THE FIRST CALL TO CHARLN. C MLIN - RETURNED CHARACTER ARRAY WITH 1 LINE. IF MLIN IS NOT C LARGE ENOUGH TO CONTAIN ALL OF THE CHARACTERS, THEY ARE C LOST. C NCML - NO. OF CHARACTERS IN THE LINE, MLIN. C BYTE MBYTA(NCBA), MCH CHARACTER*(*) MLIN IF(NCLEF.LE.0) GO TO 50 IS = NCBA - NCLEF + 1 NCMX = LEN(MLIN) J = 0 DO 40 I=IS,NCBA MCH = MBYTA(I) IF(MCH.NE.10.AND.I.LT.NCBA) GO TO 30 NCML = J NCLEF = NCBA - I GO TO 50 30 J = J + 1 IF(J.LE.NCMX) MLIN(J:J) = CHAR(MCH) 40 CONTINUE 50 RETURN END SUBROUTINE CALEND (NMONTH,NDAY,NYEAR,TAV,XDAY,FRACYR) C C DOUBLE PRECISION XN,Y,XDA DIMENSION XMOF(12) DATA XMOF(1),XMOF(2),XMOF(3),XMOF(4)/0.,31.,59.,90./ DATA XMOF(5),XMOF(6),XMOF(7),XMOF(8)/120.,151.,181.,212./ DATA XMOF(9),XMOF(10),XMOF(11),XMOF(12)/243.,273.,304.,334./ C----------------------------------------------------------------------- C THIS SUBROUTINE RETURNS THE FRACTIONAL DAY NUMBER STARTING JAN 1,1965 C (XDAY) AND FRACYR. C FRACYR=NUMBER OF YEARS IN FOUR YEAR CYCLE. C CHEAT BY REFERENCING TO 1/1/1901 C----------------------------------------------------------------------- NN=NYEAR-1 XN=NN Y=FLOAT(NN/4) XDA=365.*XN+XMOF(NMONTH)+FLOAT(NDAY)+TAV/24.+Y-1. IF((MOD(NN+1,4).EQ.0).AND.(NMONTH.GT.2)) XDA=XDA+1. FRACYR=XDA/365.25-4.*Y XDAY=XDA-23376.0D0 RETURN END SUBROUTINE EPHEM2(IHDAY, UTIM ,R0) C C INPUT DATA: C IHDAY = HOWARD DAY. C UTIM = UNIVERSAL TIME, IN HOURS. C R0 = SOLAR RADIUS IN X OR Y UNITS. C C OUTPUT DATA: C P = ANGLE BETWEEN SUN AXIS AND EARTH AXIS. C B0 = ANGLE BETWEEN SUN'S EQUATOR AND EARTH'S ORBIT. C DEC0 = DECLINATION. C SINCLN,COSCLN = SIN, COS OF ECLIPTIC LONGITUDE. C SINB0,COSB0 = SIN, COS OF B0. C SINP,COSP = SIN, COS OF P. C PHI = C FE = SCALE FACTOR = RADIANS OF ARC PER X OR Y UNIT IN THE SKY. C H0 = HOUR ANGLE. C CLONG0 = ECLIPTIC LONGITUDE. C ZD = ZENITH DISTANCE. C DOUBLE PRECISION DHDAY, DREG, DRE2 COMMON/EPHE/P,B0,DEC0,SINCLN,COSCLN,SINB0,COSB0,SINP,COSP,PHI,FE, 1 H0,CLONG0, ZD C ACCOUNT FOR LEAP SECONDS.. CONVERT U.T. TO TDT, EPHEMERIS TIME. C THIS ALGORITHM IS OFF BY MAX OF + OR - 1.3 SEC BETW. 1965 AND 1998. DREG = IHDAY DRE2 = UTIM YR = 65.D0 + DREG/365.25D0 + DRE2/8766.D0 TDTIM = DRE2 + DBLE(0.833266*YR - 17.4614)/3600.D0 C COMPUTE 365.25 DAY YEAR FROM HOWARD DAY. DHDAY = IHDAY DHDAY = DHDAY + DBLE(TDTIM)/24.D0 YR = 65.D0 + DHDAY/365.25D0 C DETERMINE THE EQUATION OF TIME COEFFICIENTS AS A FUNCTION OF YEAR. SLC = -0.137275*YR - 93.5950 S2LC = -0.869101E-2*YR + 596.946 S3LC = 0.653750E-2*YR + 3.78607 S4LC = -12.7 CLC = 0.407986E-1*YR - 432.656 C2LC = -0.3800088E-2*YR - 1.74833 C3LC = 19.3 C LONGITUDE OF ASCENDING NODE OF SOLAR EQUATOR ON THE ECLIPTIC. C ASNOLN = (0.139312E-1*YR + 74.3669)/57.29578 ASNOLN = (0.013958333*(YR + 49.997947) + 73.666667)/57.29578 C MEAN OBLIQUITY OF THE ECLIPTIC W. RESPECT TO THE MEAN C EQUATOR OF DATE. C OBEQ = (-0.129317E-3*YR + 23.452229)/57.29578 CALL GEOMLN(CLM, OBEQ, IHDAY, TDTIM) C CLM = GEOMLN(IHDAY, TDTIM) CLM2 = CLM + CLM CLM3 = CLM + CLM2 CLM4 = CLM2 + CLM2 EOT = SLC*SIN(CLM) + S2LC*SIN(CLM2) + S3LC*SIN(CLM3) + 1 S4LC*SIN(CLM4) + CLC*COS(CLM) + C2LC*COS(CLM2) + 2 C3LC*COS(CLM3) C CONVERT EQUATION OF TIME FROM SECONDS OF TIME TO RADIANS OF ANGLE. EOT = EOT*0.000072722052 C -5.2021284 IS -118.06 DEG LONGITUDE OF MT. WILSON CONVERTED C TO HOURS + 12 HOURS, THEN CONVERTED TO RADIANS. TMLO = UTIM IF(TMLO.LT.8.) TMLO = TMLO + 24. H0 = 0.26179939*TMLO - 5.2021284 + EOT RASUN = CLM - EOT COSRA = COS(RASUN) CLONG0 = ATAN((SIN(RASUN)/COSRA)/COS(OBEQ)) IF (COSRA . LT .0.) CLONG0 = CLONG0 + 3.1415927 IF (CLONG0 . LT .0.) CLONG0 = CLONG0 + 6.2831853 FE = 4.65362E-03*(1. + .016723*COS(CLONG0 + 1.35489))/R0 SINCLN = SIN(CLONG0) COSCLN = COS(CLONG0) DEC0 = ASIN(SINCLN*SIN(OBEQ)) CLOMEG = CLONG0 - ASNOLN C .12619897 IS SIN(7.25) .12721607 IS TAN(7.25) 7.25 DEG IS MAX B0. B0 = ASIN(.12619897*SIN(CLOMEG)) P = ATAN(-TAN(OBEQ)*COSCLN) + ATAN(-.12721607*COS(CLOMEG)) C BMTW IS THE LATITUDE OF MT. WILSON. BMTW = 34.216667/57.29578 COSZ = SIN(BMTW)*SIN(DEC0) + COS(BMTW)*COS(DEC0)*COS(H0) ZD = ACOS(COSZ) SINB0 = SIN(B0) COSB0 = COS(B0) SINP = SIN(P) COSP = COS(P) PHI = -ATAN(TAN(OBEQ)*COSCLN) RETURN END SUBROUTINE GEOMLN(GEOML, OBEQ, IHDAY, TDTIM) C C RETURNS THE GEOMETRIC MEAN LONGITUDE, GEOML AND THE MEAN C OBLIQUITY OF THE ECLIPTIC, OBEQ. INPUT IS HOWARD DAY, C IHDAY AND TDTIM, TERRESTRIAL DYNAMICAL TIME (U.T. + C LEAP SECONDS). C DOUBLE PRECISION DHDAY, DANT, BIGT DHDAY = IHDAY DHDAY = DHDAY + DBLE(TDTIM)/24.D0 C 23741.5 IS THE NO. OF JULIAN DAYS ELAPSED FROM 1/0.5/1900 TO C 1/1.0/1965, HOWARD DAY 0.0. DHDAY = DHDAY + 23741.5D0 C BIGT IS THE NO. OF JULIAN CENTURIES ELAPSED FROM 1/0.5/1900. BIGT = DHDAY/36525.D0 DANT = 279.696677777778D0 + 36000.768925D0*BIGT + 1 0.0003025D0*BIGT**2 DANT = DANT - DINT(DANT/360.D0)*360.D0 GEOML = DANT/57.29577951D0 DANT = 23.452294D0 - 0.0130125D0*BIGT - 1 0.00000164D0*BIGT**2 + 0.000000503D0*BIGT**3 OBEQ = DANT/57.29577951D0 RETURN END function nfloor(a) nf = nint(a) if (a.lt.0.0) nf = nf -1 nfloor = nf end