C Solution of Poisson's equation -c△u = f in the rectangle
C [XMIN,XMAX]×[YMIN,YMAX]  by Fourier series
C 近似精度は Fourier 展開係数の項数 KOSU, および, 
C 長方形上の定積分を４点公式で計算するときの辺分割数 IIB に依存する
C ここでは IB は結果のグラフを描くときのメッシュで解の近似度とは無関係
      PROGRAM POISN1
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(IB=19,KOSU=100,IIB=19)
      PARAMETER(XMIN=0,XMAX=1,YMIN=0,YMAX=1,ZMIN=-0.5,ZMAX=1,C=1)
      PARAMETER(IXMIN=0,IXMAX=799,IYMIN=0,IYMAX=599,IXD=5,IYD=6)
      REAL*8 CF(KOSU,KOSU)
      INTEGER IGCX(IB,IB),IGCY(IB,IB)
      IGX(X,Y,Z)=IXBIAS+IBREAD*(Y-YMIN)/(YMAX-YMIN)
     +                 +ISLANT*(XMAX-X)/(XMAX-XMIN)
      IGY(X,Y,Z)=ISLANT*(XMAX-X)/(XMAX-XMIN)
     +           +IBREAD*(Z-ZMIN)/(ZMAX-ZMIN)
      HX=(XMAX-XMIN)/(IB+1)
      HY=(YMAX-YMIN)/(IB+1)
      PERSP=0.8
      IBREAD=(IYMAX-IYMIN+1)/(1.0+PERSP/1.414)
      ISLANT=IBREAD*PERSP/1.414
      IXBIAS=(IXMAX-IXMIN+1)-IBREAD-ISLANT
      WRITE(*,*)'Wait a moment while calculating '
     +          //'the Fourier Coefficients.'
      CALL FRCOEF(CF,KOSU,XMIN,XMAX,YMIN,YMAX,C,IIB)
      CALL INIT(IXMIN,IYMIN,IXMAX,IYMAX)
      CALL CLS()
      X=XMIN
      IX=IGX(X,YMIN,0.0D0)
      IY=IGY(X,YMIN,0.0D0)
      IX2=IGX(X,YMAX,0.0D0)
      IY2=IGY(X,YMAX,0.0D0)
      CALL LINE(IX,IY,IX2,IY2)
      X=X+HX
      DO 300 I=0,IB-1
         XR=(X-XMIN)/(XMAX-XMIN)
         Y=YMIN
         Z=0.0D0
         IX=IGX(X,Y,Z)
         IY=IGY(X,Y,Z)
         CALL MOVE(IX,IY)
         DO 200 J=1,IB
            Y=Y+HY
            YR=(Y-YMIN)/(YMAX-YMIN)
            Z=FFOUR(XR,YR,CF,KOSU)
            IX=IGX(X,Y,Z)
            IY=IGY(X,Y,Z)
            CALL LINETO(IX,IY)
            IGCX(I+1,J)=IX
            IGCY(I+1,J)=IY
  200    CONTINUE
         CALL LINETO(IGX(X,YMAX,0.0D0),IGY(X,YMAX,0.0D0))
         X=X+HX
  300 CONTINUE
      IX=IGX(X,YMIN,0.0D0)
      IY=IGY(X,YMIN,0.0D0)
      IX2=IGX(X,YMAX,0.0D0)
      IY2=IGY(X,YMAX,0.0D0)
      CALL LINE(IX,IY,IX2,IY2)
      Y=YMIN
      IX=IGX(XMIN,Y,0.0D0)
      IY=IGY(XMIN,Y,0.0D0)
      IX2=IGX(XMAX,Y,0.0D0)
      IY2=IGY(XMAX,Y,0.0D0)
      CALL LINE(IX,IY,IX2,IY2)
      DO 500 J=1,IB
         X=XMIN
         Y=Y+HY
         Z=0.0D0
         IX=IGX(X,Y,Z)
         IY=IGY(X,Y,Z)
         CALL MOVE(IX,IY)
         DO 400 I=0,IB-1
            CALL LINETO(IGCX(I+1,J),IGCY(I+1,J))
  400    CONTINUE
         CALL LINETO(IGX(XMAX,Y,0.0D0),IGY(XMAX,Y,0.0D0))
  500 CONTINUE
      IX=IGX(XMIN,YMAX,0.0D0)
      IY=IGY(XMIN,YMAX,0.0D0)
      IX2=IGX(XMAX,YMAX,0.0D0)
      IY2=IGY(XMAX,YMAX,0.0D0)
      CALL LINE(IX,IY,IX2,IY2)
      WRITE(*,*)'Value at the center : ',FFOUR(0.5D0,0.5D0,CF,KOSU)
      PAUSE
      END
C
C 方程式の右辺の函数のフーリエ展開の係数を計算し配列CFに格納
      SUBROUTINE FRCOEF(CF,KOSU,XMIN,XMAX,YMIN,YMAX,C,IIB)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(PI=3.1415926535897932)
      REAL*8 CF(KOSU,KOSU)
      EXTERNAL G
      CC=4/(XMAX-XMIN)/(YMAX-YMIN)/C/PI/PI
      XX=(XMAX-XMIN)**2
      YY=(YMAX-YMIN)**2
      DO 200 I=1,KOSU
         DO 100 J=1,KOSU+1-I
            CF(I,J)=SQINT(G,I,J,XMIN,XMAX,YMIN,YMAX,IIB)
            CF(I,J)=CF(I,J)*CC/(I*I/XX+J*J/YY)
            IF (DABS(CF(I,J)).LT.1.0D-8) CF(I,J)=0.0D0
  100    CONTINUE
  200 CONTINUE
      RETURN
      END
C
C 配列CFに格納されたデータを係数とするフーリエ級数の和を計算する
      FUNCTION FFOUR(X,Y,CF,KOSU)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(PI=3.1415926535897932384)
      REAL*8 CF(KOSU,KOSU)
      FFOUR=0
      DO 200 I=1,KOSU
         DO 100 J=1,KOSU+1-I
            IF (CF(I,J).NE.0.0D0) THEN
               FFOUR=FFOUR+CF(I,J)*DSIN(I*PI*X)*DSIN(J*PI*Y)
            END IF
  100    CONTINUE
  200 CONTINUE
      RETURN
      END
C
C 函数F(X,Y)の長方形[XMIN,XMAX]×[YMIN,YMAX]上の二次元定積分を
C 四点公式で計算
      FUNCTION SQINT(F,I,J,XMIN,XMAX,YMIN,YMAX,IIB)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(GP=1.732050807568877/3)
      EXTERNAL F
      HX=(XMAX-XMIN)/(IIB+1)
      HY=(YMAX-YMIN)/(IIB+1)
      GX1=HX*(1-GP)/2
      GX2=HX*(1+GP)/2
      GY1=HY*(1-GP)/2
      GY2=HY*(1+GP)/2
      SQINT=0
      X=XMIN
      DO 70 K=0,IIB
         Y=YMIN
         DO 50 L=0,IIB
            F1=F(X+GX1,Y+GY1,I,J)
            F2=F(X+GX1,Y+GY2,I,J)
            F3=F(X+GX2,Y+GY1,I,J)
            F4=F(X+GX2,Y+GY2,I,J)
            SQINT=SQINT+F1+F2+F3+F4
            Y=Y+HY
   50    CONTINUE
         X=X+HX
   70 CONTINUE
      SQINT=SQINT*HX*HY/4
      RETURN
      END
C
C 方程式の右辺の函数のフーリエ展開の係数を計算するための被積分函数
      FUNCTION G(X,Y,I,J)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(XMIN=0,XMAX=1,YMIN=0,YMAX=1,PI=3.14159265358979323,
     +   C1=PI/(XMAX-XMIN),C2=PI/(YMAX-YMIN))
      G=F(X,Y)*DSIN(C1*I*(X-XMIN))*DSIN(C2*J*(Y-YMIN))
      RETURN
      END
C
C GIVEN RIGHT-HAND SIDE
      FUNCTION F(X,Y)
      REAL*8 F,X,Y
C      IF ((DABS(X-0.5).LT.0.2D0).AND.(DABS(Y-0.5).LT.0.2D0)) THEN
C       F=5
C      ELSE
C       F=0
C      END IF
      F=5
      RETURN
      END