!PROGRAM MULTPILE_CYLINDERS

!     等水深域直立任意斷面形狀柱體群受規則波作用引起相互干涉
!     本論文原始發表於1974年「複數個不透過及び透過島堤による波の散亂」，其中不透水部分理論推導簡單，適於邊界元素法初學者參考練習。
!      http://www.xn--06q65y2zay3jryqs5a29hepz.tw/chou_html/reference/1974ijima_01.pdf

      REAL,ALLOCATABLE:: X(:),Y(:),S(:)
      REAL,ALLOCATABLE:: XN(:),YN(:),GB(:,:),G(:,:),WH(:,:)
      COMPLEX,ALLOCATABLE:: FO(:),FOB(:),FF(:),FB(:),H(:,:),HB(:,:)
      COMPLEX::RS,IM,TEMP,FF0
      REAL::KH

      OPEN(UNIT=10,FILE='10.APT',FORM='UNFORMATTED',STATUS='UNKNOWN')
      OPEN(UNIT=22,FILE='WH_CYLINDERS.OUT',
     & FORM='FORMATTED',STATUS='UNKNOWN')

      IM=(0.,1.)
      PAI=3.1415927

      WRITE(*,*)  '等深域水深 Ho(M)=？？'
!      READ(*,*)    Ho 
      WRITE(*,*)  '柱體數P=？？'
!      READ(*,*)    P
      WRITE(*,*)  '矩形斷面3個呈直線為例'
      WRITE(*,*)  '間距SPAN=？？'
!      READ(*,*)    SPAN   
      WRITE(*,*)  '長邊A/Ho=？？'
!      READ(*,*)    A
      WRITE(*,*)  '短邊B/Ho=？？'
!      READ(*,*)    B
!       
!  設定柱體群尺寸(可自行調整)
      P=3
      Ho=10
      SPAN=1*Ho
      A=2*Ho;B=0.5*Ho
      NA=20;NB=4
      SA=A/NA;SB=B/NB      
      NT=2*(NA+NB)*P    ! 總節點數 採用一定元素

      ALLOCATE(GB(NT,NT),G(NT,NT),X(NT),Y(NT),S(NT),XN(NT),YN(NT))
      ALLOCATE(FO(NT),FOB(NT),FF(NT),FB(NT),H(NT,NT),HB(NT,NT))

!  設定元素座標及線段長

!  長邊分割成NA個線分，長邊分割成NB個線分
!  各矩形斷面柱體各以上方長邊右上角為起始(逆時針)

!  第1柱體
      DO I=1,NA
      II=I
      X(II)=-A/2-SPAN-SA/2-(I-1)*SA  
      Y(II)=B/2
      S(II)=SA 
      XN(II)=0
      YN(II)=1
      END DO
      DO I=1,NB
      II=I+NA
      X(II)=-A-A/2-SPAN
      Y(II)=B/2-SB/2-(I-1)*SB
      S(II)=SB 
      XN(II)=-1
      YN(II)=0
      END DO
      DO I=1,NA
      II=I+NA+NB
      X(II)=-A-A/2-SPAN+SA/2+(I-1)*SA  
      Y(II)=-B/2
      S(II)=SA 
      XN(II)=0
      YN(II)=-1
      END DO
      DO I=1,NB
      II=I+NA+NB+NA
      X(II)=-A/2-SPAN
      Y(II)=-B/2+SB/2+(I-1)*SB 
      S(II)=SB 
      XN(II)=1
      YN(II)=0
      END DO

!  第2柱體
      DO I=1,NA
      II=I+2*(NA+NB)
      X(II)=A/2-SA/2-(I-1)*SA  
      Y(II)=B/2
      S(II)=SA 
      XN(II)=0
      YN(II)=1
      END DO
      DO I=1,NB
      II=I+NA+2*(NA+NB)
      X(II)=-A/2
      Y(II)=B/2-SB/2-(I-1)*SB
      S(II)=SB 
      XN(II)=-1
      YN(II)=0
      END DO
      DO I=1,NA
      II=I+NA+NB+2*(NA+NB)
      X(II)=-A/2+SA/2+(I-1)*SA  
      Y(II)=-B/2
      S(II)=SA 
      XN(II)=0
      YN(II)=-1
      END DO
      DO I=1,NB
      II=I+NA+NB+NA+2*(NA+NB)
      X(II)=A/2
      Y(II)=-B/2+SB/2+(I-1)*SB 
      S(II)=SB 
      XN(II)=1
      YN(II)=0
      END DO

!  第3柱體
      DO I=1,NA
      II=I+4*(NA+NB)
      X(II)=A/2+A+SPAN-SA/2-(I-1)*SA  
      Y(II)=B/2
      S(II)=SA 
      XN(II)=0
      YN(II)=1
      END DO
      DO I=1,NB
      II=I+NA+4*(NA+NB)
      X(II)=A/2+SPAN
      Y(II)=B/2-SB/2-(I-1)*SB
      S(II)=SB 
      XN(II)=-1
      YN(II)=0
      END DO
      DO I=1,NA
      II=I+NA+NB+4*(NA+NB)
      X(II)=A/2+SPAN+SA/2+(I-1)*SA  
      Y(II)=-B/2
      S(II)=SA 
      XN(II)=0
      YN(II)=-1
      END DO
      DO I=1,NB
      II=I+NA+NB+NA+4*(NA+NB)
      X(II)=A+SPAN
      Y(II)=-B/2+SB/2+(I-1)*SB 
      S(II)=SB 
      XN(II)=1
      YN(II)=0
      END DO

!  波浪資訊
        WRITE(*,*) '輸入入射波週期(秒) T=？？'
        READ(*,*) T
!        WRITE(*,*) '輸入入射波有義波高(m) H=？？'
!        READ(*,*) H
        WRITE(*,*)  '輸入入射波角度(°)'   ! 與X軸所呈角度
        READ(*,*) OO

        OI=OO/180.*PAI        !   入射波角度轉換成弧度
        SHG=4*PAI*PAI*Ho/T/T/9.8

        CALL SHG_KH(SHG,KH)   !   牛頓近似法求KH

        WRITE(*,*) '週期T=',T, '  入射波角度OO=',OO

!  柱體上入射波勢函數及法線方向導函數

        DO I=1,NT
        TEMP=CEXP(-IM*KH*(X(I)*COS(OI)+Y(I)*SIN(OI)))
        FO(I)=-IM*TEMP
        FOB(I)=-TEMP*(COS(OI)*XN(I)+SIN(OI)*YN(I))*KH
        END DO

!  2維邊界積分方程式

      CALL G2_HANKEL(GB,G,M,N,X,Y,S,XN,YN,KH,H,HB)  !結果存於[H] 

!  柱體斷面各節點散射波速度勢

      DO I=1,NT
      FB(I)=-FOB(I)
      END DO

      DO I=1,NT
      TEMP=0
      DO J=1,NT
      TEMP=TEMP+H(I,J)*FOB(J)
      END DO
      FF(I)=-TEMP    
      END DO

!  柱體群週邊水面波高分佈
!  設定範圍 DISxDIS(可自行調整)
!  各柱體內事後必須去除

      DIS=5*Ho
      XL=DIS*(2*PAI*Ho/KH)
      YL=XL
      NX=50;NY=50

      ALLOCATE(WH(NX,NY))

      DO IX=1,2*NX
      XX=-XL+XL/NX*(IX-1)
      DO IY=1,2*NY
      YY=YL-YL/NY*(IY-1)

      FF0=-IM*KH*CEXP(-IM*KH*(XX*COS(OI)+YY*SIN(OI)))

      TEMP=0
      DO J=1,NT
      R=SQRT((XX-X(J))**2+(YY-Y(J))**2)
      RK=KH*R
      CALL BESJNS(RK,0,BJO,ILL)
      IF(ILL.NE.0) WRITE(*,1) ILL
1     FORMAT(1X,'ILL G2 =',I7)
      CALL BESJNS(RK,1,BJ1,ILL)
      IF(ILL.NE.0) WRITE(*,1) ILL
      CALL BESYNS(RK,0,BYO,ILL)
      IF(ILL.NE.0) WRITE(*,1) ILL
      CALL BESYNS(RK,1,BY1,ILL)
      IF(ILL.NE.0) WRITE(*,1) ILL
      RN=((X(J)-XX)*YN(J)-(Y(J)-YY)*XN(J))/R
      RS=0.5*IM*KH*(BJ1+IM*BY1)*S(J)*RN
      RT=-0.5*IM*(BJO+IM*BYO)*S(J)
      TEMP=TEMP+RT*FB(J)-RS*FF(J)
      END DO

      WH(IX,IY)=CABS(FF0+TEMP)    !水面波高分佈

      END DO   ! IY
      END DO   ! IX

      WRITE(22,3) ((WH(IX,IY),IY=1,NY),IX=1,NX)
3     FORMAT(50F5.2)

      STOP  
      END

!  副程式

      SUBROUTINE G2_HANKEL(GB,G,M,N,X,Y,S,XN,YN,KH,H,HB)

! 2維平面邊界積分方程式
! 結果存於H ，H為複數陣列
! KI：=1 順時針積分；=2 逆時針積分

      REAL::GB(M,M),G(M,M),X(N),Y(N),S(N),XN(N),YN(N),KH
      COMPLEX::H(M,M),HB(M,M),IM,RS,RT

      IM=(0.,1.)
      PAI=3.1415927

      DO I=1,N
      DO J=1,N

      IF(I.NE.J) THEN
      R=SQRT((X(I)-X(J))**2+(Y(I)-Y(J))**2)
      RK=KH*R

      CALL BESJNS(RK,0,BJO,ILL)
      IF(ILL.NE.0) WRITE(*,1) ILL
1     FORMAT(1X,'ILL G2 =',I7)
      CALL BESJNS(RK,1,BJ1,ILL)
      IF(ILL.NE.0) WRITE(*,1) ILL
      CALL BESYNS(RK,0,BYO,ILL)
      IF(ILL.NE.0) WRITE(*,1) ILL
      CALL BESYNS(RK,1,BY1,ILL)
      IF(ILL.NE.0) WRITE(*,1) ILL

      RN=((X(J)-X(I))*YN(J)-(Y(J)-Y(I))*XN(J))/R
      RS=0.5*IM*KH*(BJ1+IM*BY1)*S(J)*RN
      RT=-0.5*IM*(BJO+IM*BYO)*S(J)
      END IF
 
      KI=2        !預設

      IF(I.EQ.J) THEN
      IF(KI.EQ.1) RS=(-1.,0.) ! KI：=1 順時針積分     
      IF(KI.EQ.2) RS=(1.,0.)  ! KI：=2 逆時針積分    
      RT=(ALOG(KH*S(J)/4.)-0.42278-IM*PAI/2.)/PAI*S(J)
      END IF

      GB(I,J)=REAL(RS)
      G(I,J)=AIMAG(RS)

      HB(I,J)=RT

      END DO
      END DO

      REWIND 15
      WRITE(15) ((GB(I,J),J=1,N),I=1,N)
      WRITE(15) ((G(I,J),J=1,N),I=1,N)

      CALL CMINVS(GB,G,M,N)

      DO I=1,N
      DO J=1,N
      RS=(0.,0.)
      DO K=1,N
      RS=RS+(GB(I,K)+IM*G(I,K))*HB(K,J)
      END DO
      H(I,J)=RS     !  [H]^[G}
      END DO
      END DO

      RETURN
      END

      SUBROUTINE CMINVS(A,B,M,N)

!   解複數連立方程式(應用逆矩陣)A+iB
!   原數據依實數部虛數部序，存檔於15
!   結果：實數矩陣存於A，實數矩陣存於B
!   14：作業檔
 
      REAL::A(M,M),B(M,M)

      IM=(0.,1.)
      EPS=1E-7     !有關pivot精度，視必要可適確調整
      REWIND 15
      READ(15) ((A(I,J),J=1,N),I=1,N)
      CALL MINVS(A,M,N,EPS,ILL)
      IF(ILL.NE.0) WRITE(*,10) ILL
10    FORMAT(1X,'ILL CS.1 ='I7)
      READ(15) ((B(I,J),J=1,N),I=1,N)
      CALL SEKI(A,B,M,N,N,N)
      REWIND 14
      WRITE(14)  ((A(I,J),J=1,N),I=1,N)
      CALL SEKI(B,A,M,N,N,N)
      REWIND 15
      READ(15) ((A(I,J),J=1,N),I=1,N)

      DO I=1,N
      DO J=1,N
      A(I,J)=A(I,J)+B(I,J)
      END DO
      END DO

      CALL MINVS(A,M,N,EPS,ILL)
      IF(ILL.NE.0) WRITE(*,11) ILL
11    FORMAT(1X,'ILL CS.2 =',I7)
      REWIND 14
      READ(14) ((B(I,J),J=1,N),I=1,N)
      CALL SEKI(B,A,M,N,N,N)

      DO I=1,N
      DO J=1,N
      B(I,J)=-B(I,J)
      END DO
      END DO

      RETURN
      END

      SUBROUTINE SHG_KH(SHG,HK)

! **利用牛頓近似法，已知無因次週頻率shg從分散關係式求 kh**

	IF((SHG-10.).LE.0) GOTO 2
	XY=SHG
	GO TO 6
2	IF((SHG-1.0).GE.0) GOTO 4
	X=SQRT(SHG)
	GO TO 5
4	X=SHG
5	COTHX=1./TANH(X)
	XY=X-(X-SHG*COTHX)/(1.+SHG*(COTHX**2-1.))
	E=1.-XY/X
	X=XY      
	IF((ABS(E)-0.0005).GE.0) GOTO 5
6	hk=XY
        RETURN
        END

      SUBROUTINE BESJNS(X,N,BJN,ILL)

      BJN=0.0
      IF(N.GE.0) GOTO 20
      ILL=1
      RETURN
20    IF(X.LT.0) GOTO 30
      IF(X.EQ.0) GOTO 21
      IF(X.LT.0) GOTO 31
21    IF(N.GT.1) GOTO 30
      ILL=0
      IF(N.EQ.0) BJN=1.0
      RETURN
30    ILL=2
      RETURN
31    IF((X-15.0).GT.0) GOTO 34
      NTEST=20.0+10.0*X-X*X/3.0
      GOTO 36
34    NTEST=110.0+X/2.0
36    IF((N-NTEST).LT.0) GOTO 40
38    ILL=3
      RETURN
40    ILL=0
!
! -- *** COMPUTE STARTING VALUE OF M *** --
!
      N1=0
      IF(X.LT.1.0)   GOTO 54
      IF(X.LT.10.0)  GOTO 53
      IF(X.LT.50.0)  GOTO 52
      IF(X.LT.100.0) GOTO 51
      M=X+35.0
      GOTO 55
51    M=1.1*X+25.0
      GOTO 55
52    M=1.3*X+15.0
      N1=N+10
      GOTO 55
53    M=2.0*X+8.0
      N1=N+7
      GOTO 55
54    M=4.0*X+6.0
55    M=MAX0(M,N1)
!
! --- ** SET F(M),F(M-1) ** ---
!
      FM1=1.0E-30
      FM=0.0
      ALPHA=0.0
      IF((M-(M/2)*2).NE.0) GOTO 120
      JT=-1
      GOTO 130
120   JT=1
130   M2=M-2
      DO K=1,M2
      MK=M-K
      FMK=FLOAT(MK)
      BMK=2.0*FMK*FM1/X-FM
      FM=FM1
      FM1=BMK
      IF((MK-N-1).NE.0) GOTO 150
      BJN=BMK
150   JT=-JT
      S=1+JT
      ALPHA=ALPHA+BMK*S
      END DO
      BMK=2.0*FM1/X-FM
      IF(N.NE.0) GOTO 180
       BJN=BMK
180   ALPHA=ALPHA+BMK
      BJN=BJN/ALPHA
      RETURN
      END

      SUBROUTINE BESYNS(X,N,BYN,ILL)

      DOUBLE PRECISION TERM
      BYN=-1.0E30
      IF(N.LT.0) GOTO 180
      ILL=0
      IF(X.LT.0) GOTO 190
      IF(X.EQ.0) GOTO 11
      IF(X.GT.0) GOTO 20
11    IF(N.GT.1) GOTO 190
      RETURN
20    PI=3.141592653
      IF((X-4.0).LE.0) GOTO 40
      T=4.0/X
      T=T*T
      P0=((((-0.0000037043*T+0.0000173565)*T-0.0000487613)*T
     &  +0.0001734300)*T-0.0017530620)*T+0.3989422793

      Q0=((((0.0000032312*T-0.0000142078)*T+0.0000342468)*T
     &  -0.0000869791)*T+0.0004564324)*T-0.0124669441

      P1=((((0.0000042414*T-0.0000200920)*T+0.0000580759)*T
     &  -0.0002232030)*T+0.0029218256)*T+0.3989422819

      Q1=((((-0.0000036594*T+0.0000162200)*T-0.0000398708)*T
     &  +0.0001064741)*T-0.0006390400)*T+0.0374008364

      A=SQRT(2.0*PI)
      B=4.0*A
      P0=A*P0
      Q0=B*Q0/X
      P1=A*P1
      Q1=B*Q1/X
      A=X-PI/4.0
      B=SQRT(2.0/(PI*X))
      Y0=B*(P0*SIN(A)+Q0*COS(A))
      Y1=B*(-P1*COS(A)+Q1*SIN(A))
      GOTO 90
!
! --- ** COMPUTE Y0,Y1  DO X <= 4 ** ---
!
40    XX=X/2.0
      X2=XX*XX
      T=ALOG(XX)+0.5772156649
      SUM=0.0
      TERM=T
      Y0=T

      DO L=1,15
      IF((L-1).EQ.0) GOTO 60
      SUM=SUM+1.0/FLOAT(L-1)
60    FLL=L
      TS=T-SUM
      TERM=(TERM*(-X2)/(FLL*FLL))*(1.0-1.0/(FLL*TS))
      Y0=Y0+TERM
      END DO
      TERM=XX*(T-0.5)

      SUM=0.0
      Y1=TERM
      DO L=2,16
      SUM=SUM+1.0/FLOAT(L-1)
      FLL=L
      FL1=FLL-1.0
      TS=T-SUM
      TERM=(TERM*(-X2)/(FL1*FLL))*((TS-0.5/FLL)/(TS+0.5/FL1))
      Y1=Y1+TERM
      END DO


      PI2=2.0/PI
      Y0=PI2*Y0
      Y1=-PI2/X+PI2*Y1
!
!只求Y0 或 Y1

90    IF((N-1).GT.0) GOTO 130

      IF(N.EQ.0) GOTO 120
      BYN=Y1
      RETURN
120   BYN=Y0
      RETURN   !只求Y0 或 Y1

!                 求 YN

130   YA=Y0
      YB=Y1
      K=1
140   T=FLOAT(2*K)/X
      YC=T*YB-YA
      K=K+1
      IF((K-N).EQ.0) GOTO 160
      YA=YB
      YB=YC
      GOTO 140
160   BYN=YC
      RETURN
180   ILL=1
      RETURN
190   ILL=2
      RETURN
      END

      SUBROUTINE SEKI(A,B,M,NA,NK,NB)

!  矩陣相乘
!  [A],[B]宣告(m, m)
!  [A]=[A][B]
!  (na,nb)=(na,nk)x(nk,nb)
!  [A]資料被取代

      REAL:: A(M,M),B(M,M)
      REAL:: C(NK)
      DO I=1,NA
      DO K=1,NK
      C(K)=A(I,K)
      END DO
      DO J=1,NB
      R=0
      DO K=1,NK
      R=R+C(K)*B(K,J)
      END DO
      A(I,J)=R
      END DO
      END DO

      WRITE(*,*) 'SEKI'
      RETURN
      END

      SUBROUTINE MINVS(A,KO,N,EPS,ILL)

      REAL::A(KO,KO)
      INTEGER:: X(KO)
      LOGICAL B
      ILL=0
      IF((KO.GE.N).AND.(N.GE.2).AND.(EPS.GT.0.0)) GOTO 1
      ILL=30000
      RETURN
1     DO I=1,N
      X(I)=I
      END DO
      DO K=1,N
      M=K
      W=0.0
      DO I=K,N
      ABSA=ABS(A(I,K))
      IF(ABSA.LE.W) GOTO 20
      W=ABSA
      M=I
20    END DO
      IF(M.EQ.K) GOTO 50
      L=X(M)
      X(M)=X(K)
      X(K)=L
      DO J=1,N
      W=A(K,J)
      A(K,J)=A(M,J)
      A(M,J)=W
      END DO
50    IF(ABS(A(K,K)).GE.EPS) GOTO 60
      ILL=K
      RETURN
60    P=1.0/A(K,K)
      DO  J=1,N
      A(K,J)=A(K,J)*P
      END DO
      DO I=1,N
      T=-A(I,K)
      B=(I.NE.K).AND.(T.NE.0.0)
      IF(.NOT.B) GOTO 100
      DO J=1,N
      A(I,J)=A(I,J)+A(K,J)*T
      END DO
100   A(I,K)=P*T
      END DO
      A(K,K)=P
      END DO
      DO I=1,N
120   IF(X(I).EQ.I) GOTO 140
      K=X(I)
      DO J=1,N
      W=A(J,I)
      A(J,I)=A(J,K)
      A(J,K)=W
      END DO
      L=X(I)
      X(I)=X(K)
      X(K)=L
      GOTO 120
140   END DO

      WRITE(*,*) 'MINVS'
      RETURN
      END
!
        SUBROUTINE WA(A,B,NA,N,NS)

!      [A]N,NS + [B]N,NS =====> [B]N,NS

        REAL::A(NA,NA),B(NA,NA)
        DO I=1,N
        DO J=1,NS
        B(I,J)=A(I,J)+B(I,J)
        END DO
        END DO
        WRITE(*,*) 'WA'
        RETURN
        END
!
        SUBROUTINE SA(A,B,NA,N,NS)

!      [A]N,NS - [B]N,NS =====> [B]N,NS

        REAL::A(NA,NA),B(NA,NA)
        DO I=1,N
        DO J=1,NS
        B(I,J)=A(I,J)-B(I,J)
        END DO
        END DO
        WRITE(*,*) 'SA'
        RETURN
        END