SUBROUTINE PSOLVE(Z,PSI,WORK,NX,NY,DX) ! ---------------------------------------------------------------------- ! Poisson Solver ! * Periodic boundary conditions in x ! (final 2 columns equivalent to first two columns) ! * Fixed boundary values in y (Z=0) ! * NX must be (2**n1)+2 ! * NY must be (2**n2)+1 ! * WORK must have dimension at least 2*(NX-2)*(NY-1) ! ! Calls 'Numerical Recipes' N-Dimensional Fast Fourier Transform ! Routine FOURN ! ! Based on the 'Numerical Recipes' 1D real transform routine ! REALFFT. ! ! Data is extended to double the channel width to allow for ! components with wavelengths which are half integer not ! integer multiples of the width BUT compacted to half the ! channel length by utilising both the Real and Imaginary ! parts of the input array. ! ---------------------------------------------------------------------- ! ! !--DECLARE PARAMETERS REAL PI, PI2 PARAMETER( PI=3.1415927,PI2=6.2831854 ) ! !---DECLARE VARIABLES PASSED THROUGH ARGUMENT LIST INTEGER NX, NY REAL Z(NX,NY),PSI(NX,NY),WORK(*) REAL DX ! !---DECLARE LOCAL VARIABLES INTEGER IX,IX1,IX2,IY,IY1,IY2,KOUNT1,KOUNT2,KOUNTR,KOUNTY, & NCENX,NCENY,NELEM,NUMX,NUMY REAL AK1,AK2,AL,F1R,F1I,F2R,F2I,FSR,FDR,FDI,FSI,FACTX,FACTY, & FCTR1,FCTR2,H1R,H1I,H2R,H2I,HSR,HDR,HSI,HDI,CS,SN,PSITOP, & PSIBOT,SCALEF,THETA INTEGER NN(2) ! NUMX=NX-2 NUMY=NY-1 NELEM=2*NUMX*NUMY ! NN(1)=NUMX/2 NN(2)=NUMY*2 ! ! ---------------------------------------------------------------------- ! Fill work array WORK with data from Z ! KOUNTR=1 ! ! Force Z=0 on lower boundary ! DO IX=1,NUMX WORK(KOUNTR)=0 KOUNTR=KOUNTR+1 ENDDO ! DO IY=2,NUMY DO IX=1,NUMX WORK(KOUNTR)=+Z(IX,IY) KOUNTR=KOUNTR+1 ENDDO ENDDO ! ! Force Z=0 on upper boundary ! DO IX=1,NUMX WORK(KOUNTR)=0 KOUNTR=KOUNTR+1 ENDDO ! ! Copy 'Reflected' Z to WORK ! DO IY=NUMY,2,-1 DO IX=1,NUMX WORK(KOUNTR)=-Z(IX,IY) KOUNTR=KOUNTR+1 ENDDO ENDDO ! ! ---------------------------------------------------------------------- ! Perform Fourier Transform ! CALL FOURN(WORK,NN,2,1) ! ! ---------------------------------------------------------------------- ! Divide fourier coefficients of Z by -(total-wavenumber)**2 ! to give coefficients of PSI ! ! NOTE: The coefficient of PSI for (total-wavenumber)=0 is ! set to 0. ! ! x-wavenumber (k) = AKi = IXi*FACTX (i=1 or 2) ! y-wavenumber (l) = AL = IY *FACTY ! ! This program works on TWO frequencies at a time and ! only calculates a little over half the updated ! coefficients. The remainder are filled by symmetry. ! FACTX=PI2/(DX*REAL(NUMX)) FACTY=PI2/(DX*REAL(NN(2))) THETA=PI2/REAL(NUMX) NCENX=NN(1)/2 NCENY=NN(2)/2 ! DO IX1=0,NCENX IX2=NN(1)-IX1 CS=COS(IX1*THETA) SN=SIN(IX1*THETA) AK1=IX1*FACTX AK2=IX2*FACTX ! DO IY=0,NCENY AL=IY*FACTY KOUNTY=IY*NUMX+1 KOUNT1=KOUNTY+2*IX1 KOUNT2=KOUNTY+2*IX2 IF (IX1.EQ.0) KOUNT2=KOUNT1 ! ! NOTE: H2 is minus the value found in WORK. ! H1R=WORK(KOUNT1) H1I=WORK(KOUNT1+1) H2R=-WORK(KOUNT2) H2I=-WORK(KOUNT2+1) ! ! Calculate coefficients corresponding to true frequencies ! using compacted complex coefficients H1 and H2 ! HSR=H1R+H2R HDR=H1R-H2R HSI=H1I+H2I HDI=H1I-H2I ! F1R=0.5*(HSR+CS*HSI+SN*HDR) F1I=0.5*(HDI-CS*HDR+SN*HSI) F2R=0.5*(HSR-CS*HSI-SN*HDR) F2I=0.5*(-HDI-CS*HDR+SN*HSI) ! ! Calculate scaling factors -1/(total wavenumber)**2 ! IF (AK1.NE.0.OR.AL.NE.0) THEN FCTR1=-1/(AK1*AK1+AL*AL) ELSE FCTR1=0. ENDIF ! FCTR2=-1/(AK2*AK2+AL*AL) ! F1R=F1R*FCTR1 F1I=F1I*FCTR1 F2R=F2R*FCTR2 F2I=F2I*FCTR2 ! FSR=F1R+F2R FDR=F1R-F2R FSI=F1I+F2I FDI=F1I-F2I ! ! Place modified H1 and H2 back in WORK array. ! WORK(KOUNT1)=0.5*(FSR-CS*FSI+SN*FDR) WORK(KOUNT1+1)=0.5*(FDI+CS*FDR+SN*FSI) WORK(KOUNT2)=-0.5*(FSR+CS*FSI-SN*FDR) WORK(KOUNT2+1)=-0.5*(-FDI+CS*FDR+SN*FSI) ! ENDDO ENDDO ! ! Fill remainder of WORK array using symmetry ! DO IY1=1,NCENY-1 IY2=NN(2)-IY1 KOUNT1=IY1*NUMX KOUNT2=IY2*NUMX DO IX=1,NUMX KOUNT1=KOUNT1+1 KOUNT2=KOUNT2+1 WORK(KOUNT2)=-WORK(KOUNT1) ENDDO ENDDO ! ! ---------------------------------------------------------------------- ! Inverse Fourier Transform ! CALL FOURN(WORK,NN,2,-1) ! ! ---------------------------------------------------------------------- ! Extract PSI from WORK ! ! NOTE: When performing an inverse transform , produces ! values which are a factor of NN(1)*NN(2) too large ! KOUNTR=NUMX+1 SCALEF=1./REAL((NN(1)*NN(2))) DO IY=2,NUMY DO IX=1,NUMX PSI(IX,IY)=WORK(KOUNTR)*SCALEF KOUNTR=KOUNTR+1 ENDDO ENDDO ! ! ---------------------------------------------------------------------- ! Apply boundary conditions ! ! Add linear term in y based on boundary values of PSI ! DO IX=1,NUMX PSITOP=PSI(IX,NY) PSIBOT=PSI(IX,1) SCALEF=(PSITOP-PSIBOT)/REAL(NUMY) DO IY=2,NUMY PSI(IX,IY)=PSIBOT+SCALEF*REAL(IY-1)+PSI(IX,IY) ENDDO ENDDO ! ! Apply periodicity in x ! DO IY=1,NY PSI(NX-1,IY)=PSI(1,IY) PSI(NX,IY)=PSI(2,IY) ENDDO ! RETURN END ! ! ====================================================================== ! SUBROUTINE FOURN(DATA,NN,NDIM,ISIGN) ! ! Taken from Numerical Recipes. ! ! DATA is replace by its NDIM-ensional discrete Fourier transform, if ISIGN ! is input as 1. NN is an integer array of length NDIM, containing the ! lengths of each dimension (number of complex values), which MUST all be ! powers of 2. DATA is a real array of length twice the product of these ! lengths, in which the data are stored as in a multidimensional complex ! FORTRAN array. If ISIGN is input as -1, DATA is replaced by its inverse ! transform times the product of the lengths of all dimensions. ! !---DECLARE VARIABLES PASSED THROUGH ARGUMENT LIST INTEGER ISIGN,NDIM INTEGER NN(NDIM) REAL DATA(*) ! !---DECLARED DOUBLE PRECISION FOR TRIGONOMETRIC RECURRENCES DOUBLE PRECISION WR,WI,WPR,WPI,WTEMP,THETA ! !---DECLARE REMAINING LOCAL VARIABLES INTEGER I1,I2,I2REV,I3,I3REV,IBIT,IDIM,IFP1,IFP2,IP1,IP2,IP3,N, & NPREV,NREM,NTOT,K1,K2 REAL TEMPR,TEMPI NTOT=1 ! !---COMPUTE TOTAL NUMBER OF COMPLEX VALUES DO IDIM=1,NDIM NTOT=NTOT*NN(IDIM) ENDDO NPREV=1 ! !---MAIN LOOP OVER THE DIMENSIONS DO IDIM=1,NDIM N=NN(IDIM) NREM=NTOT/(N*NPREV) IP1=2*NPREV IP2=IP1*N IP3=IP2*NREM I2REV=1 ! !---THIS IS THE BIT REVERSAL SECTION OF THE ROUTINE DO I2 = 1,IP2,IP1 IF (I2.LT.I2REV) THEN DO I1 = I2,I2+IP1-2,2 DO I3 = I1,IP3,IP2 I3REV=I2REV+I3-I2 TEMPR=DATA(I3) TEMPI=DATA(I3+1) DATA(I3)=DATA(I3REV) DATA(I3+1)=DATA(I3REV+1) DATA(I3REV)=TEMPR DATA(I3REV+1)=TEMPI ENDDO ENDDO ENDIF IBIT=IP2/2 DO WHILE ((IBIT.GE.IP1).AND.(I2REV.GT.IBIT)) I2REV=I2REV-IBIT IBIT=IBIT/2 ENDDO I2REV=I2REV+IBIT ENDDO ! !---DANIELSON-LANCZOS SECTION OF THE ROUTINE STARTS HERE IFP1=IP1 DO WHILE (IFP1.LT.IP2) IFP2=2*IFP1 ! !---INITIALISE FOR THE TRIG. RECURRENCE THETA=ISIGN*6.28318530717959D0/(IFP2/IP1) WPR=-2.D0*DSIN(0.5D0*THETA)**2 WPI=DSIN(THETA) WR=1.D0 WI=0.D0 DO I3=1,IFP1,IP1 DO I1=I3,I3+IP1-2,2 DO I2=I1,IP3,IFP2 ! !---DANIELSON-LANCZOS FORMULA K1=I2 K2=K1+IFP1 TEMPR=SNGL(WR)*DATA(K2)-SNGL(WI)*DATA(K2+1) TEMPI=SNGL(WR)*DATA(K2+1)+SNGL(WI)*DATA(K2) DATA(K2)=DATA(K1)-TEMPR DATA(K2+1)=DATA(K1+1)-TEMPI DATA(K1)=DATA(K1)+TEMPR DATA(K1+1)=DATA(K1+1)+TEMPI ENDDO ENDDO ! !---TRIGONOMETRIC RECURRENCE WTEMP=WR WR=WR*WPR-WI*WPI+WR WI=WI*WPR+WTEMP*WPI+WI ENDDO IFP1=IFP2 ENDDO NPREV=N*NPREV ENDDO RETURN END