program MOC C******************************************************************************** C c Multiple Object Convolution C C This routine is designed to take an object file C convolve with each of a set of the input files, i.e. the convolution c images to produce a set of corresponding images. C C Written: Julian Christou C DATE: September 30, 1996 C C******************************************************************************** c c Dimensioning Of Arrays c c******************************************************************************** parameter(mm=16384) real*4 data(mm), obj(mm), psf(mm) integer get_input character dfile*80, ofile*80, pfile*80, errmsg*80 integer lud, luo, lup, dtype, dtypeo data lud/55/, luo/56/, lup/57/ c******************************************************************************** c c Logfile Information c c******************************************************************************** open(unit=77,file='moc.log',form='formatted') write(77,10000) 10000 format(6x,'MOC Log File',/) c******************************************************************************* c c Open Data Files (NOTE ALL ARE IN IRAF FORMAT) c c******************************************************************************* jer=0 write(6,100) 100 format(3x,'Multiple Object Convolution Routine',/,/, * 6x,'Enter PSF Data Filename: ',$) 101 istat=get_input(dfile,' ',nlen) call imfile_open(dfile,lud,ier,nx,ny,nz,dtype) if(ier.ne.0) then jer=jer + 1 if(jer.eq.5) stop write(6,105) jer 105 format(6x,i3,': Open Error - Enter New Filename: ',$) goto 101 endif write(77,106) dfile 106 format(/,3x,'Data File: ',a40) write(6,110) nx, ny, nz write(77,110) nx, ny, nz 110 format(/,6x,'File is ',i3,' x ',i3,' x ',i4,' Pixels') write(6,120) dtype write(77,120) dtype 120 format(6x,'Data Type = ',i5,/) if(nx.eq.ny) then nd = nx ndd = nd * nd nd2 = nd / 2 nd21 = nd2 + 1 else write(6,130) 130 format(6x,'Size Error - Data Frames are not Square') stop endif c******************************************************************************* c c Open Object File (NOTE ALL ARE IN IRAF FORMAT) c c******************************************************************************* jer=0 write(6,200) 200 format(6x,'Enter Object Data Filename: ',$) 201 istat=get_input(ofile,' ',nlen) call imfile_open(ofile,luo,ier,nxo,nyo,nzo,dtypeo) if(ier.ne.0) then jer=jer + 1 if(jer.eq.5) stop write(6,105) jer goto 201 endif write(77,206) ofile 206 format(/,3x,'Object File: ',a40) write(6,110) nxo, nyo, nzo write(77,110) nxo, nyo, nzo write(6,120) dtypeo write(77,120) dtypeo if(nxo.eq.nyo) then ndo = nxo nddo = ndo * ndo nd2o = ndo / 2 nd21o = nd2o + 1 else write(6,130) stop endif call imfile_read_2d(luo, obj, nxo, nyo) C do j = 1, nddo C print *, j, obj(j) C enddo C****************************************************************************** c c Enter Output Filename Parameters c c***************************************************************************** write(6,300) 300 format(6x,'Enter Output Filename: ',$) istat = get_input(pfile,' ',nlen) write(77,310) pfile 310 format(/,6x,'Ouput File : ',a20) write(6,301) 301 format(6x,'Enter Beginning Frame Number: ',$) read(5,*) n1 write(6,302) 302 format(6x,'Enter Ending Frame Number: ',$) read(5,*) n2 write(77,303) n1, n2 write(6,303) n1,n2 303 format(6x,'Processing Frames',i3,' to ',i3) nz = n2 - n1 + 1 call iraf_write_open_3d(pfile,lup,nd,nd,nz,dtype) c***************************************************************************** c c Processing Loop c c***************************************************************************** do ifr = n1, n2 call imfile_read_3d(lud, data, nx, ny, ifr) C print *, ifr C do j = 1, nddo C print *, ifr, j, psf(j) C enddo c**************************************************************************** c c Convolve Data Frame By Object c c**************************************************************************** call con(data, obj, psf, nd, nd) call iraf_write_3d(psf, nd, nd, ifr, lup) enddo c****************************************************************************** c c TERMINATE ROUTINE & Close Output Arry c c****************************************************************************** call imclos(lup,ier) if(ier.ne.0) then call imemsg(ier,errmsg) write(6,500) ier, errmsg 500 format(/,6x,i3,':',a80,/) endif stop end c****************************************************************************** c c NECESSARY SUBROUTINES & FUNCTIONS c c****************************************************************************** c c get_input c c***************************************************************************** function get_input(result,prompt,i) character*(*) result,prompt integer i integer get_input print 30,prompt 30 format (' ',(a),$) read (*,100) result 100 format((a)) i = len(result) do while ((result(i:i) .eq. ' ') .and. (i.ne.0)) i = i-1 end do get_input = 0 return end c**************************************************************** c c imfile_open c c**************************************************************** subroutine imfile_open(file,im,ierr,nx,ny,nz,dtype) C THIS SUBROUTINE IS SET UP TO OPEN THE IRAF FILE WHICH CONTAINS THE C IMAGES OF THE 2-D IR DATA SETS. C JULIAN CHRISTOU C 881128 C MODIFIED TO HANDLE DON'S 3D DATA SETS C JULIAN CHRISTOU C 890921 integer ier,axlen(7),naxis,dtype character file*80,errmsg*80 C OPEN IRAF IMAGE ierr=0 call imopen(file,1,im,ier) if(ier.ne.0) then call imemsg(ier,errmsg) write(6,110) errmsg 110 format(/,6x,a80,/) ierr=1 goto 500 endif C READ IN THE IMAGE HEADER 501 call imgsiz(im,axlen,naxis,dtype,ier) if(ier.ne.0) then call imemsg(ier,errmsg) write(6,110) errmsg goto 501 endif nx = axlen(1) ny = axlen(2) nz = axlen(3) 500 return end c******************************************************************************* c c imfile_read_2d c c**************************************************************** subroutine imfile_read_2d(im,dummy,nx,ny) C SUBROUTINE TO READ A SINGLE IMAGE FROM AN IRAF FILE OPENED WITH C imfile_open. real*4 dummy(nx*ny) character errmsg*80 mx = nx my = ny C print *, nx,ny jer = 0 call imgs2r(im,dummy,1,mx,1,my,ier) if(ier.ne.0) then jer = jer + 1 if(jer.eq.5) stop call imemsg(ier,errmsg) write(6,100) errmsg 100 format(/,6x,a80,/) stop endif do j = 8192, 8292 print *, j, dummy(j) enddo return end c***************************************************************** c c imfile_read_3d c c**************************************************************** subroutine imfile_read_3d(im,dummy,nx,ny,ifile) C SUBROUTINE TO READ THE ifile FRAME FROM THE IRAF FILE OPENED C BY imfile_open. real*4 dummy(nx*ny) character errmsg*80 jer = 0 call imgs3r(im,dummy,1,nx,1,ny,ifile,ifile,ier) if(ier.ne.0) then jer = jer + 1 if(jer.eq.5) stop call imemsg(ier,errmsg) write(6,100) errmsg 100 format(/,6x,a80,/) stop endif do j = 8192, 8292 print *, j, dummy(j) enddo return end c***************************************************************** c c con c c**************************************************************** subroutine con(data,obj,psf,nx,ny) c This subroutine computes the Fourier transforms of the input data array c and the object array and takes the Fourier quotien, inverts and returns c the ouput to the main routine real*4 data(nx*ny), dummy(16384), obj(nx*ny), odummy(16384) real*4 psf(16384) integer*4 ndim(2) complex*8 cdummy(8320), cdum(16384), ocdummy(8320) complex*8 ocdum(16384) equivalence (dummy,cdummy),(odummy,ocdummy) nd1 = nx * ny print *, nx, ny, nd1 nd21 = nd /2 +1 ndim(1) = nd1 ndim(2) = nd2 c Copy Data into working array do i = 1, nd1 dummy(i) = data(i) enddo c Compute Data FFT call reass(dummy,nd,nd) call four2(dummy,ndim,2,-1,0) c Pad out Data to full transform call chtcf(cdummy,cdum,nd,nd,1) call creass(cdum,nd,nd) c Copy Obj into working array do i = 1, nd1 odummy(i) = obj(i) enddo c Compute Obj FFT call reass(odummy,nd,nd) call four2(odummy,ndim,2,-1,0) c Pad out Obj to full transform call chtcf(ocdummy,ocdum,nd,nd,1) call creass(ocdum,nd,nd) c Compute Quotient of Complex arrays do i = 1, nd1 cdum(i) = cdum(i) * ocdum(i) enddo c Replace into half-plane call creass(cdum,nd,nd) call chtcf(cdummy,cdum,nd,nd,-1) c Inverse Transform to new image call four2(dummy,ndim,2,1,-1) call reass(dummy,nd,nd) c Copy working array back into data array do i = 1, nd1 psf(i) = dummy(i) / float(nd1) enddo do j = 8192, 8292 print *, j, psf(j) enddo return end c***************************************************************** c c chtcf c c**************************************************************** Subroutine chtcf(cdata1,cdata2,mx,my,idir) C Routine Complex Half plane To Complex Full plane. This takes a C complex half plane image CDATA1 as outputted from FOUR2 with iform=0 C and creates a full complex image CDATA2. C N.B. It is assumed the DC point is at the top LEFT corner! C IDIR: +1 Complex/2[cdata1] to Complex[cdata2] C -1 Complex[cdata2] to Complex/2[cdata1] C Julian C Christou 850301 complex cdata1(mx/2+1,my),cdata2(mx,my) if(idir.eq.1) then do 10 i=1,my do 10 j=1,mx/2+1 10 cdata2(j,i)=cdata1(j,i) do 20 j=1,mx/2-1 cdata2(mx+1-j,1)=conjg(cdata2(j+1,1)) 20 cdata2(mx+1-j,my/2+1)=conjg(cdata2(j+1,my/2+1)) do 30 i=1,my/2-1 do 30 j=1,mx/2-1 cdata2(mx+1-j,my+1-i)=conjg(cdata2(j+1,i+1)) 30 cdata2(mx+1-j,i+1)=conjg(cdata2(j+1,my+1-i)) return elseif(idir.eq.-1) then do 40 i=1,my do 40 j=1,mx/2+1 40 cdata1(j,i)=cdata2(j,i) return endif write(10,100) idir 100 format(/,'********** INVALID DIRECTION FLAG IN CHTCF',i5, * ' **********',/) return end c***************************************************************** c c creass c c**************************************************************** Subroutine creass(data,mx,my) C Routine to REASSemble a complex 2-D array as outputed from an FFT to place C the DC point (or equivalent) at the center. Clockwise quadrants 1234 C become 3412. DATA is a REAL array of size MX x MY. C Julian C Christou 850301 C modified by Ron. 11apr85 complex*8 data(mx,my),a mx2=mx/2 my2=my/2 do 10 i=1,my2 do 10 j=1,mx2 a=data(j,i) data(j,i)=data(j+mx2,i+my2) data(j+mx2,i+my2)=a a=data(j+mx2,i) data(j+mx2,i)=data(j,i+my2) data(j,i+my2)=a 10 continue return end c***************************************************************** c c reass c c**************************************************************** Subroutine reass(data,mx,my) C Routine to REASSemble a 2-D array as outputed from an FFT to place C the DC point (or equivalent) at the center. Clockwise quadrants 1234 C become 3412. DATA is a REAL array of size MX x MY. C Julian C Christou 850301 C modified by Ron. 11apr85 real*4 data(mx,my),a mx2=mx/2 my2=my/2 do 10 i=1,my2 do 10 j=1,mx2 a=data(j,i) data(j,i)=data(j+mx2,i+my2) data(j+mx2,i+my2)=a a=data(j+mx2,i) data(j+mx2,i)=data(j,i+my2) data(j,i+my2)=a 10 continue return end c***************************************************************** c c four2 c c**************************************************************** SUBROUTINE FOUR2 (DAT,N,NDIM,ISIGNE,IFORM) C COOLEY-TUKEY FAST FOURIER TRANSFORM IN ANSI BASIC FORTRAN. C MULTI-DIMENSIONAL TRANSFORM, EACH DIMENSION A POWER OF TWO, C COMPLEX OR REAL DATA. C TRANSFORM(K1,K2,...) = SUM(DAT(J1,J2,...)*EXP(ISIGNE*2*PI*SQRT(-1) C *((J1-1)*(K1-1)/N(1)+(J2-1)*(K2-1)/N(2)+...))), SUMMED FOR ALL C J1 AND K1 FROM 1 TO N(1), J2 AND K2 FROM 1 TO N(2), ETC. FOR ALL C NDIM SUBSCRIPTS. NDIM MUST BE POSITIVE AND ALL N(KDIM) POWERS OF C TWO. ISIGNE IS + OR -1. LET NTOT = N(1)*N(2)*...*N(NDIM). THEN A C -1 TRANSFORM FOLLOWED BY A +1 ONE (OR VICE VERSA) RETURNS NTOT C TIMES THE ORIGINAL DATA. IFORM = 1, 0 OR -1, AS DATA IS C COMPLEX, REAL OR THE FIRST HALF OF A COMPLEX ARRAY. TRANSFORM C VALUES ARE RETURNED TO ARRAY DAT. THEY ARE COMPLEX, REAL OR C THE FIRST HALF OF A COMPLEX ARRAY, AS IFORM = 1, -1 OR 0. C THE TRANSFORM OF A REAL ARRAY (IFORM = 0) DIMENSIONED N(1) BY N(2) C BY ... WILL BE RETURNED IN THE SAME ARRAY, NOW CONSIDERED TO C BE COMPLEX OF DIMENSIONS N(1)/2+1 BY N(2) BY .... NOTE THAT IF C IFORM = 0 OR -1, N(1) MUST BE EVEN, AND ENOUGH ROOM MUST BE C RESERVED. THE MISSING VALUES MAY BE OBTAINED BY COMPLEX CONJUGA- C TION. THE REVERSE TRANSFORMATION, OF A HALF COMPLEX ARRAY DIMEN- C SIONED N(1)/2+1 BY N(2) BY ..., IS ACCOMPLISHED BY SETTING IFORM C TO -1. IN THE N ARRAY, N(1) MUST BE THE TRUE N(1), NOT N(1)/2+1. C THE TRANSFORM WILL BE REAL AND RETURNED TO THE INPUT ARRAY. C RUNNING TIME IS PROPORTIONAL TO NTOT*LOG2(NTOT), RATHER THAN C THE NAIVE NTOT**2. C WRITTEN BY NORMAN BRENNER, MIT, JUNE 1968. SEE-- IEEE AUDIO C TRANSACTIONS, JUNE 1967 AND JUNE 1969, SPECIAL FFT ISSUES. DIMENSION DAT(524288), N(ndim) NTOT=1 DO 10 KDIM=1,NDIM IF (N(KDIM)) 110,110,10 10 NTOT=NTOT*N(KDIM) IF (IFORM) 70,20,20 20 NREM=NTOT DO 60 KDIM=1,NDIM NREM=NREM/N(KDIM) NPREV=NTOT/(N(KDIM)*NREM) NCURR=N(KDIM) IF (KDIM-1+IFORM) 30,30,40 30 NCURR=NCURR/2 40 CALL BITRV (DAT,NPREV,NCURR,NREM) CALL DFT2 (DAT,NPREV,NCURR,NREM,ISIGNE) IF (KDIM-1+IFORM) 50,50,60 50 CALL FIXRL (DAT,N(1),NREM,ISIGNE,IFORM) NTOT=(NTOT/N(1))*(N(1)/2+1) 60 CONTINUE RETURN 70 NTOT=(NTOT/N(1))*(N(1)/2+1) NREM=1 DO 100 JDIM=1,NDIM KDIM=NDIM+1-JDIM NCURR=N(KDIM) IF (KDIM-1) 80,80,90 80 NCURR=NCURR/2 CALL FIXRL (DAT,N(1),NREM,ISIGNE,IFORM) NTOT=NTOT/(N(1)/2+1)*N(1) 90 NPREV=NTOT/(N(KDIM)*NREM) CALL BITRV (DAT,NPREV,NCURR,NREM) CALL DFT2 (DAT,NPREV,NCURR,NREM,ISIGNE) 100 NREM=NREM*N(KDIM) RETURN 110 WRITE (6,120) KDIM,N(KDIM) 120 FORMAT (20H0ERROR IN FOUR2. N(,I2,4H) = ,I5,45H IS ZERO OR NEGATI *VE. NO TRANSFORM WAS DONE.) RETURN END c***************************************************************** c c fixrl c c**************************************************************** SUBROUTINE FIXRL (DAT,N,NREM,ISIGNE,IFORM) C FOR IFORM = 0, CONVERT THE TRANSFORM OF A DOUBLED-UP REAL ARRAY, C CONSIDERED COMPLEX, INTO ITS TRUE TRANSFORM. SUPPLY ONLY THE C FIRST HALF OF THE COMPLEX TRANSFORM, AS THE SECOND HALF HAS C CONJUGATE SYMMETRY. FOR IFORM = -1, CONVERT THE FIRST HALF C OF THE TRUE TRANSFORM INTO THE TRANSFORM OF A DOUBLED-UP REAL C ARRAY. N MUST BE EVEN. C USING COMPLEX NOTATION AND SUBSCRIPTS STARTING AT ZERO, THE C TRANSFORMATION IS-- C DIMENSION DAT(N/2,NREM) C ZSTP = EXP(ISIGNE*2*PI*I/N) C DO 10 I2=0,NREM-1 C DAT(0,I2) = CONJ(DAT(0,I2))*(1+I) C DO 10 I1=1,N/4 C Z = (1+(2*IFORM+1)*I*ZSTP**I1)/2 C I1CNJ = N/2-I1 C DIF = DAT(I1,I2)-CONJ(DAT(I1CNJ,I2)) C TEMP = Z*DIF C DAT(I1,I2) = (DAT(I1,I2)-TEMP)*(1-IFORM) C 10 DAT(I1CNJ,I2) = (DAT(I1CNJ,I2)+CONJ(TEMP))*(1-IFORM) C IF I1=I1CNJ, THE CALCULATION FOR THAT VALUE COLLAPSES INTO C A SIMPLE CONJUGATION OF DAT(I1,I2). C REMOVE THE NEXT TWO CARDS IF NO DOUBLE PRECISION FEATURE EXISTS. C THEY SERVE TO INCREASE ACCURACY OF THE SINE AND COSINE VALUES. C DOUBLE PRECISION TWOPI,THETA,SIN,SINTH,ZSTPR,ZSTPI,ZR,ZI,TEMP C SIN(THETA)=DSIN(THETA) DIMENSION DAT(524288) TWOPI=6.283185*FLOAT(ISIGNE) IP0=2 IP1=IP0*(N/2) IP2=IP1*NREM IF (IFORM) 10,70,70 C PACK THE REAL INPUT VALUES (TWO PER COLUMN) 10 J1=IP1+1 DAT(2)=DAT(J1) IF (NREM-1) 70,70,20 20 J1=J1+IP0 I2MIN=IP1+1 DO 60 I2=I2MIN,IP2,IP1 DAT(I2)=DAT(J1) J1=J1+IP0 IF (N-2) 50,50,30 30 I1MIN=I2+IP0 I1MAX=I2+IP1-IP0 DO 40 I1=I1MIN,I1MAX,IP0 DAT(I1)=DAT(J1) DAT(I1+1)=DAT(J1+1) 40 J1=J1+IP0 50 DAT(I2+1)=DAT(J1) 60 J1=J1+IP0 70 DO 80 I2=1,IP2,IP1 TEMPR=DAT(I2) DAT(I2)=DAT(I2)+DAT(I2+1) 80 DAT(I2+1)=TEMPR-DAT(I2+1) IF (N-2) 200,200,90 90 THETA=TWOPI/FLOAT(N) SINTH=SIN(THETA/2.) ZSTPR=-2.*SINTH*SINTH ZSTPI=SIN(THETA) ZR=(1.-ZSTPI)/2. ZI=(1.+ZSTPR)/2. IF (IFORM) 100,110,110 100 ZR=1.-ZR ZI=-ZI 110 I1MIN=IP0+1 I1MAX=IP0*(N/4)+1 DO 190 I1=I1MIN,I1MAX,IP0 DO 180 I2=I1,IP2,IP1 I2CNJ=I2+IP0*(N/2)-2*(I1-1) IF (I2-I2CNJ) 150,120,120 120 IF (ISIGNE*(2*IFORM+1)) 130,140,140 130 DAT(I2+1)=-DAT(I2+1) 140 IF (IFORM) 170,180,180 150 DIFR=DAT(I2)-DAT(I2CNJ) DIFI=DAT(I2+1)+DAT(I2CNJ+1) TEMPR=DIFR*ZR-DIFI*ZI TEMPI=DIFR*ZI+DIFI*ZR DAT(I2)=DAT(I2)-TEMPR DAT(I2+1)=DAT(I2+1)-TEMPI DAT(I2CNJ)=DAT(I2CNJ)+TEMPR DAT(I2CNJ+1)=DAT(I2CNJ+1)-TEMPI IF (IFORM) 160,180,180 160 DAT(I2CNJ)=DAT(I2CNJ)+DAT(I2CNJ) DAT(I2CNJ+1)=DAT(I2CNJ+1)+DAT(I2CNJ+1) 170 DAT(I2)=DAT(I2)+DAT(I2) DAT(I2+1)=DAT(I2+1)+DAT(I2+1) 180 CONTINUE TEMP=ZR-.5 ZR=ZSTPR*TEMP-ZSTPI*ZI+ZR 190 ZI=ZSTPR*ZI+ZSTPI*TEMP+ZI 200 IF (IFORM) 270,210,210 C UNPACK THE REAL TRANSFORM VALUES (TWO PER COLUMN) 210 I2=IP2+1 I1=I2 J1=IP0*(N/2+1)*NREM+1 GO TO 250 220 DAT(J1)=DAT(I1) DAT(J1+1)=DAT(I1+1) I1=I1-IP0 J1=J1-IP0 230 IF (I2-I1) 220,240,240 240 DAT(J1)=DAT(I1) DAT(J1+1)=0. 250 I2=I2-IP1 J1=J1-IP0 DAT(J1)=DAT(I2+1) DAT(J1+1)=0. I1=I1-IP0 J1=J1-IP0 IF (I2-1) 260,260,230 260 DAT(2)=0. 270 RETURN END c***************************************************************** c c bitrv c c**************************************************************** SUBROUTINE BITRV (DAT,NPREV,N,NREM) C PERMUTE THE DATA BY ^BIT REVERSAL^. C DIMENSION DAT(NPREV,N,NREM) C COMPLEX DAT C EXCHANGE DAT(J1,J4REV,J5) WITH DAT(J1,J4,J5) FOR ALL J1 FROM 1 C TO NPREV, ALL J4 FROM 1 TO N (WHICH MUST BE A POWER OF TWO), AND C ALL J5 FROM 1 TO NREM. J4REV-1 IS THE BIT REVERSAL OF J4-1. E.G. C SUPPOSE N = 32. THEN FOR J4-1 = 10011, J4REV-1 = 11001, ETC. DIMENSION DAT(524288) IP0=2 IP1=IP0*NPREV IP4=IP1*N IP5=IP4*NREM I4REV=1 C I4REV = 1+(J4REV-1)*IP1 DO 60 I4=1,IP4,IP1 C I4 = 1+(J4-1)*IP1 IF (I4-I4REV) 10,30,30 10 I1MAX=I4+IP1-IP0 DO 20 I1=I4,I1MAX,IP0 C I1 = 1+(J1-1)*IP0+(J4-1)*IP1 DO 20 I5=I1,IP5,IP4 C I5 = 1+(J1-1)*IP0+(J4-1)*IP1+(J5-1)*IP4 I5REV=I4REV+I5-I4 C I5REV = 1+(J1-1)*IP0+(J4REV-1)*IP1+(J5-1)*IP4 TEMPR=DAT(I5) TEMPI=DAT(I5+1) DAT(I5)=DAT(I5REV) DAT(I5+1)=DAT(I5REV+1) DAT(I5REV)=TEMPR 20 DAT(I5REV+1)=TEMPI C ADD ONE TO THE HIGH ORDER BIT OF J4REV-1 WITH DOWNWARD CARRY. 30 IP2=IP4/2 40 IF (I4REV-IP2) 60,60,50 50 I4REV=I4REV-IP2 IP2=IP2/2 IF (IP2-IP1) 60,40,40 60 I4REV=I4REV+IP2 RETURN END c***************************************************************** c c dft2 c c**************************************************************** SUBROUTINE DFT2 (DAT,NPREV,N,NREM,ISIGNE) C DISCRETE FOURIER TRANSFORM OF LENGTH N. IN-PLACE COOLEY-TUKEY C ALGORITHM, BIT-REVERSED TO NORMAL ORDER, SANDE-TUKEY PHASE SHIFTS. C DIMENSION DAT(NPREV,N,NREM) C DAT(J1,K4,J5) = SUM(DAT(J1,J4,J5)*EXP(ISIGNE*2*PI*I*(J4-1)* C (K4-1)/N)), SUMMED OVER J4 = 1 TO N FOR ALL J1 FROM 1 TO NPREV, C K4 FROM 1 TO N AND J5 FROM 1 TO NREM. N MUST BE A POWER OF TWO. C METHOD--LET IPREV TAKE THE VALUES 1, 2 OR 4, 4 OR 8, ..., N/16, C N/4, N. THE CHOICE BETWEEN 2 OR 4, ETC., DEPENDS ON WHETHER N IS C A POWER OF FOUR. DEFINE IFACT = 2 OR 4, THE NEXT FACTOR THAT C IPREV MUST TAKE, AND IREM = N/(IFACT*IPREV). THEN-- C DIMENSION DAT(NPREV,IPREV,IFACT,IREM,NREM) C COMPLEX DAT C DAT(J1,J2,J3,J4,J5) = SUM(DAT(J1,J2,J3,J4,J5)*EXP(ISIGNE*2*PI*I* C (K3-1)*((J3-1)/IFACT+(J2-1)/(IFACT*IPREV)))), SUMMED OVER J3 = 1 C TO IFACT FOR ALL J1 FROM 1 TO NPREV, J2 FROM 1 TO IPREV, K3 FROM C 1 TO IFACT, J4 FROM 1 TO IREM AND J5 FROM 1 TO NREM. THIS IS C A PHASE-SHIFTED DISCRETE FOURIER TRANSFORM OF LENGTH IFACT. C FACTORING N BY FOURS SAVES ABOUT TWENTY FIVE PERCENT OVER FACTOR- C ING BY TWOS. DATA MUST BE BIT-REVERSED INITIALLY. DIMENSION DAT(524288) C REMOVE THE NEXT TWO CARDS IF NO DOUBLE PRECISION FEATURE EXISTS. C THEY SERVE TO INCREASE ACCURACY OF THE SINE AND COSINE VALUES. C DOUBLE PRECISION SIN,TWOPI,THETA,SINTH,WSTPR,WSTPI,WR,WI,TEMP C SIN(THETA)=DSIN(THETA) TWOPI=6.283185*FLOAT(ISIGNE) IP0=2 IP1=IP0*NPREV IP4=IP1*N IP5=IP4*NREM IP2=IP1 NPART=N 10 IF (NPART-2) 70,30,20 20 NPART=NPART/4 GO TO 10 C DO A FOURIER TRANSFORM OF LENGTH TWO 30 IF (IP2-IP4) 40,170,170 40 IP3=IP2*2 C IP2=IP1*IPROD C IP3=IP2*IFACT DO 50 I1=1,IP1,IP0 C I1 = 1+(J1-1)*IP0 DO 50 I5=I1,IP5,IP3 C I5 = 1+(J1-1)*IP0+(J4-1)*IP3+(J5-1)*IP4 I3A=I5 I3B=I3A+IP2 C I3 = 1+(J1-1)*IP0+(J2-1)*IP1+(J3-1)*IP2+(J4-1)*IP3+(J5-1)*IP4 TEMPR=DAT(I3B) TEMPI=DAT(I3B+1) DAT(I3B)=DAT(I3A)-TEMPR DAT(I3B+1)=DAT(I3A+1)-TEMPI DAT(I3A)=DAT(I3A)+TEMPR 50 DAT(I3A+1)=DAT(I3A+1)+TEMPI C DO A FOURIER TRANSFORM OF LENGTH FOUR (FROM BIT REVERSED ORDER) 60 IP2=IP3 70 IF (IP2-IP4) 80,170,170 80 IP3=IP2*4 THETA=TWOPI/FLOAT(IP3/IP1) SINTH=SIN(THETA/2.) WSTPR=-2.*SINTH*SINTH C COS(THETA)-1, FOR ACCURACY. WSTPI=SIN(THETA) WR=1. WI=0. DO 160 I2=1,IP2,IP1 C I2 = 1+(J2-1)*IP1 IF (I2-1) 100,100,90 90 W2R=WR*WR-WI*WI W2I=2.*WR*WI W3R=W2R*WR-W2I*WI W3I=W2R*WI+W2I*WR 100 I1MAX=I2+IP1-IP0 DO 150 I1=I2,I1MAX,IP0 C I1 = 1+(J1-1)*IP0 DO 150 I5=I1,IP5,IP3 C I5 = 1+(J1-1)*IP0+(J4-1)*IP3+(J5-1)*IP4 I3A=I5 I3B=I3A+IP2 I3C=I3B+IP2 I3D=I3C+IP2 C I3 = 1+(J1-1)*IP0+(J2-1)*IP1+(J3-1)*IP2+(J4-1)*IP3+(J5-1)*IP4 IF (I2-1) 120,120,110 C APPLY THE PHASE SHIFT FACTORS 110 TEMPR=DAT(I3B) DAT(I3B)=W2R*TEMPR-W2I*DAT(I3B+1) DAT(I3B+1)=W2R*DAT(I3B+1)+W2I*TEMPR TEMPR=DAT(I3C) DAT(I3C)=WR*TEMPR-WI*DAT(I3C+1) DAT(I3C+1)=WR*DAT(I3C+1)+WI*TEMPR TEMPR=DAT(I3D) DAT(I3D)=W3R*TEMPR-W3I*DAT(I3D+1) DAT(I3D+1)=W3R*DAT(I3D+1)+W3I*TEMPR 120 T0R=DAT(I3A)+DAT(I3B) T0I=DAT(I3A+1)+DAT(I3B+1) T1R=DAT(I3A)-DAT(I3B) T1I=DAT(I3A+1)-DAT(I3B+1) T2R=DAT(I3C)+DAT(I3D) T2I=DAT(I3C+1)+DAT(I3D+1) T3R=DAT(I3C)-DAT(I3D) T3I=DAT(I3C+1)-DAT(I3D+1) DAT(I3A)=T0R+T2R DAT(I3A+1)=T0I+T2I DAT(I3C)=T0R-T2R DAT(I3C+1)=T0I-T2I IF (ISIGNE) 130,130,140 130 T3R=-T3R T3I=-T3I 140 DAT(I3B)=T1R-T3I DAT(I3B+1)=T1I+T3R DAT(I3D)=T1R+T3I 150 DAT(I3D+1)=T1I-T3R TEMP=WR WR=WSTPR*WR-WSTPI*WI+WR 160 WI=WSTPR*WI+WSTPI*TEMP+WI GO TO 60 170 RETURN END c------------------------------------------------------------------------------------- c c iraf_write_open_3d c c------------------------------------------------------------------------------------- subroutine iraf_write_open_3d(file,im,nx,ny,nz,dtype) c c Routine which creates and opens a 3D IRAF file for writing into c integer naxis, dtype, axlen(7), ier, naxis character errmsg*80, file*80 C CREATE IRAF IMAGE axlen(1) = nx axlen(2) = ny axlen(3) = nz naxis=3 dtype=6 call imcrea( file, axlen, naxis, dtype, ier) if(ier.ne.0) then call imemsg(ier,errmsg) write(6,10100) errmsg 10100 format(/,6x,a80,/) endif C OPEN IRAF IMAGE FOR WRITING call imopen( file, 3, im, ier) if(ier.ne.0) then call imemsg(ier,errmsg) write(6,10100) errmsg endif return end c***************************************************************************** c c iraf_write_3d c c**************************************************************************** subroutine iraf_write_3d(data,nx,ny,iz,im) c c Routine to write out a single frame to a 3d IRAF data file which has c been already opened. c real*4 data(nx*ny) character errmsg*80 c print*, im, nx, ny, iz C WRITE IMAGE TO IRAF IMAGE SECTION call imps3r( im, data, 1, nx, 1, ny, iz, iz, ier) if(ier.ne.0) then call imemsg(ier,errmsg) write(6,10100) ier, errmsg 10100 format(/,6x,i3,':',a80,/) endif return end