C**********************************************************************
C*  Estimation of Roll Damping for conventional cargo ship & barge    *
C*                                                                    *
C*  SA_RDPM.BAS                                                       *
C*  since 2007.09, developed & coded by Toru KATAYAMA, Dr.Eng         *
C*                                      Osaka Prefecture University   *
C**********************************************************************
C
      PROGRAM MAIN
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION NUE,L,KAI,L0,L1,L2,NABLA,KW,KAI1,LTHETA,MGA,KN
     &    ,LO,LR,KG,LB1,LB2,LB3,LB2S0,M1,M2,M3,M4,M5,M6,M7,M8,LOMEGA
      INTEGER BWOSM,TYPE_BK,RDP,CSDMIN,SSBK1,SSBK2,P
C     SC:number of SS, PT:number of data in a SS, NNRA:number of wave length for culculation
C     NNFN:number of Fn for culculation, NNKAI:number of wave direction for culculation
C     NNRD:number of roll angles for culculation
      PARAMETER SC=200,PT=200,NNRA=70,NNFN=20,NNKAI=35,NNRD=30
C
      DIMENSION  CSD(SC),NUMPT(SC),HT(0:SC,0:PT),HB(0:SC,0:PT)
      DIMENSION  THETA1(NNRD),RAM(NNRA),FUN(NNFN),KAI(NNKAI)
      DIMENSION  XX(SC,PT*2),YY(SC,PT*2)
      DIMENSION  X(SC),HO(SC),SIG(SC),BX(SC),DX(SC)
      DIMENSION  PERI(NNRA,NNFN,NNKAI),BE2(SC)
      DIMENSION  RMAX1(2),V(2),CR(SC),CR1(SC)
      DIMENSION  X1(30),YI(30),XBK(0:SC),RATIO(SC)
      DIMENSION  X11(SC),Y11(SC),WX(SC),WY(SC)
      DIMENSION  BKT(SC),BKB(SC),NUBK(SC),BBKHAT(SC)
      DIMENSION  BBYD(SC),SGM(SC),BB0(SC),TT0(SC),WLHB(SC),WLAR(SC)
      DIMENSION  GLEN(SC)
C
C
C
C******
C     if BWOSM=0 then wave making component is not culculated
      BWOSM=0
C     kinematic viscous coefficient
      NUE=1.14E-06
C     circular constant
      PAI=3.14159
C     Acceleration of gravity
      G=9.8
      RO=102.0
C
C
       OPEN(1,FILE='OSM2.RDP')
C****** data of hull form
C      Lpp, Breadth, Draft
       READ(1,*)L,B,D
C      number of square stations,  the number of mid-ship
       READ(1,*)NUMCS,IMID
       DO I=1 , NUMCS
C        SS
         READ(1,*)CSD(I)
C        number of datum in a square station (ss)
         READ(1,*)NUMPT(I)
         DO J=1 , NUMPT(I)
C          hight from Base Line, breadth from Center Line(please input from BL)
           READ(1,*)HT(I,J),HB(I,J)
         END DO
       END DO
C****** data for bilgekeel
C      breadth of bilgekeel, afr-end(ss), forward-end(ss)  (in Japan, AP=ss0, FP=ss10)
       READ(1,*)BBK,XBK1,XBK2
C      Selection of the calculation method (1: simplified, 2:integration)
       READ(1,*)TYPE_BK
       IF( TYPE_BK.NE.1)then
C        the aft-end and the forward-end of the ss number of the section where the bilge keel exists on hull's ss.
         READ(1,*)SSBK1,SSBK2
         DO I=SSBK1 , SSBK2
C          The attached transverse position of the bilge-keel from BL
           READ(1,*)BKB(I)
         END DO
       END IF
C****** eddy making component for two type ships (1: conventuonal ship, 2:barge)
       READ(1,*)RDP
C****** OGD=(Draft-KG)/Draft
       READ(1,*)OGD
C****** number of roll amplitude for culculation, roll amplitude(degree)
       READ(1,*)NTHETA
       DO I=1 , NTHETA
         READ(1,*)THETA1(I)
       END DO
C****** Fn for culculation
       READ(1,*)NFN
       DO I=1 , NFN
         READ(1,*)FUN(I)
       END DO
C****** wavelength/Lpp for culculation
       READ(1,*)NRA
       DO I=1 , NRA
         READ(1,*)RAM(I)
       END DO
C****** wave direction for culculation (degree)
       READ(1,*)NKAI
       DO I=1 , NKAI
         READ(1,*)KAI(I)
       END DO
       CLOSE(1)
C===== end of loading the data file
C
C=====in order to culculate bilge keel component by integration method, following culculation is required
        DO I=1 , NUMCS
          HT(I,NUMPT(I)+1)=999
          HB(I,NUMPT(I)+1)=HB(I,NUMPT(I))
        END DO
        DO I=1 , NUMCS
          DO J=1 , NUMPT(I)
              XX(I,NUMPT(I)-J+1)= HB(I,J)
            XX(I,NUMPT(I)+J-1)=-HB(I,J)
            YY(I,NUMPT(I)-J+1)= HT(I,J)
            YY(I,NUMPT(I)+J-1)= HT(I,J)
          END DO
        END DO
C
C====== culculation of data for hull form  =====
C======half breadth, water line area, and girth length
        WLHT=D
        DO I=1 , NUMCS
          SUM=0.0
 10    CONTINUE
          IF( HT(I,1).GT.WLHT)THEN
            WLHB(I)=0.
            WLAR(I)=0.
            GLEN(I)=0.
            GOTO 10
          END IF
          IF( HT(I,1).EQ.WLHT)THEN
            WLHB(I)=HB(I,1)
            WLAR(I)=0.
            GLEN(I)=HB(I,1)
          END IF
          DO J=1 , NUMPT(I)
 11    CONTINUE
            IF( (HT(I,J).LT.WLHT) .AND. (WLHT.LE.HT(I,J+1)))THEN
              GOTO 12
            END IF
            SUM=SUM+(HT(I,J+1)-HT(I,J))*(HB(I,J+1)+HB(I,J))/2
            GOTO 14
C      ====== calculation on each water line
 12    CONTINUE
          AA=(HB(I,J+1)-HB(I,J))/(HT(I,J+1)-HT(I,J))
          BB=(-HB(I,J+1)*HT(I,J)+HB(I,J)*HT(I,J+1))/(HT(I,J+1)-HT(I,J))
            WLHB(I)=AA*WLHT+BB
            WLAR(I)=SUM+(WLHT-HT(I,J))*(WLHB(I)+HB(I,J))/2
            GOTO 15
 14    CONTINUE
          END DO
 15    CONTINUE
        END DO
C      ===== Calculation of a section modulus
        DO I=1 , NUMCS
          IF( HT(1,1)  .LT.D)THEN
            CSDMIN=1
          END IF
          IF( (HT(I-1,1) .GT.D) .AND. (HT(I,1).GE.D))THEN
            GOTO21
          END IF
          IF( (HT(I-1,1).GE.D) .AND. (HT(I,1).LT.D))THEN
            CSDMIN=I
          END IF
          IF( (HT(I,1).GE.D) .AND. (HT(I-1,1).LT.D))THEN
            NUMH0=I-1
            GOTO 22
          END IF
          HB2=HB(I,1)
          DO J=2 , NUMPT(I)
            IF(HT(I,J).GT.WLHT)THEN
              GOTO 18
            END IF
            IF(HB(I,J).GT.HB2)THEN
              HB2=HB(I,J)
            END IF
          END DO
 18    CONTINUE
          IF( WLHB(I).GT.HB2)THEN
            HB2=WLHB(I)
          END IF
          IF( WLHB(I).EQ.0)THEN
            BBYD(I)=HB2/(WLHT-HT(I,1))
            GOTO 19
          END IF
            BBYD(I)=WLHB(I)/(WLHT-HT(I,1))
 19    CONTINUE
          IF( WLHB(I).EQ.0)THEN
            SGM(I)=WLAR(I)/(HB2*(WLHT-HT(I,1)))
            GOTO 20
          END IF
            SGM(I)=WLAR(I)/(WLHB(I)*(WLHT-HT(I,1)))
 20    CONTINUE
            TT0(I)=(WLHT-HT(I,1))/D
          IF( TT0(I).LT.0)THEN
            TT0(I)=0
          END IF
          IF( WLHB(I).EQ.0)THEN
            BB0(I)=HB2/B*2
            GOTO 21
          END IF
            BB0(I)=WLHB(I)/B*2
 21    CONTINUE
        END DO
        NUMH0=NUMCS
 22    CONTINUE
C
       N=NUMCS-CSDMIN+1
       DO I=1 , N
         X(I)  =CSD(CSDMIN-1+I)
         HO(I) =BBYD(CSDMIN-1+I)
         SIG(I)=SGM(CSDMIN-1+I)
         BX(I) =BB0(CSDMIN-1+I)*B
         DX(I) =TT0(CSDMIN-1+I)*D
       END DO
C
C====== culculation of sidplacement
       ODD=0.0
         if (mod(NUMCS,2).LT.0.5) then
          ODD=1.0
       end if
         DO J=1 , NUMCS-2-ODD , 2
         Q0=CSD(J)
         Q1=CSD(J+1)
         Q2=CSD(J+2)
         QQ=(Q0+Q2)/2
         HH=(Q2-Q0)/2*L/10/3
         L0=4*(QQ-Q1)*(QQ-Q2)/(Q0-Q1)/(Q0-Q2)+1
         L1=4*(QQ-Q0)*(QQ-Q2)/(Q1-Q0)/(Q1-Q2)
         L2=4*(QQ-Q0)*(QQ-Q1)/(Q2-Q0)/(Q2-Q1)+1
         SUM3=SUM3+(L0*WLAR(J)+L1*WLAR(J+1)+L2*WLAR(J+2))*HH
       END DO
       IF( ODD.LT.0.5)THEN
         GOTO 1404
       END IF
       J=NUMCS-1
       HH=(CSD(J+1)-CSD(J))/2*L/10
       SUM3=SUM3+(WLAR(J)+WLAR(J+1))*HH
 1404  NABLA=2*SUM3
C
C====== Calculation of a block coefficient and an inside crossection coefficient
       CB=NABLA/L/B/D
       CM=2*WLAR(IMID)/B/D
C
C======Calculation of the encounter period which calculates a roll damping (an incident-wave period, Fn, angle of wave direction)
       DO K=1 , NKAI
         DO J=1 , NFN
           DO I=1 , NRA
             OMG=SQRT(2*PAI*G/RAM(I)/L)
             KW=2*PAI/RAM(I)/L
             V2=FUN(J)*SQRT(L*G)
             KAI1=KAI(K)*PAI/180
             MGA=OMG-KW*V2*COS(KAI1)
             MGA=ABS(MGA)
             PERI(I,J,K)=2*PAI/MGA
           END DO
         END DO
       END DO
C
C===== Calculation of the attachment height of a bilge keel, and an index of a data number
       IF( TYPE_BK.NE.1)THEN
         DO I=SSBK1 , SSBK2
           J=1
 3     CONTINUE
           NUBK(I)=NUMPT(I)-J+1
           IF( HB(I,J).LE.BKB(I))THEN
             J=J+1
             GOTO 3
           END IF
           BKT(I)=(HT(I,J)-HT(I,J-1))*(BKB(I)-HB(I,J-1))
     &            /(HB(I,J)-HB(I,J-1))+HT(I,J-1)
        END DO
      END IF
C
C*********************************************************************
C
        OPEN(2,FILE='OSMB44.CSV')
        NO=0
        DO ITHETA=1 , NTHETA
          THTA=THETA1(ITHETA)*3.14159/180
          DO IKAI=1 , NKAI
            DO IFN=1 , NFN
              DO IRAM=1 , NRA
                NO=NO+1
C
                OMEGA=6.28318/PERI(IRAM,IFN,IKAI)
C
C      *********** Frictional Component ************
C      Ref. H.Kato (JZK. No.102) and  S.Tamiya  et al (JZK. No.132)
C
                SF=L*(1.7*D+CB*B)
                RF=((0.887+0.145*CB)*(1.7*D+CB*B)-2*OGD*D)/3.1415
                BFHAT=0.787*SF*RF**2.*SQRT(OMEGA*NUE*B/19.6133)
     &     /(NABLA*B**2.)*(1.+4.1*FUN(IFN)/OMEGA*SQRT(9.80665/L))
C
C
C     **********  Wave Making Component **********
C     Ref.   Y.Ikeda  et al  (JZK. No.143)
C
         GUZAID=OMEGA**2.*D/9.80665
         A1=1+GUZAID**(-1.2)*EXP(-2*GUZAID)
         A2=0.5+GUZAID**(-1)*EXP(-2*GUZAID)
         LOMEGA=OMEGA*FUN(IFN)*SQRT(L/9.80665)
         SIT=20.*(LOMEGA-0.3)

         TAH=(EXP(SIT)-EXP(-SIT))/(EXP(SIT)+EXP(-SIT))

         BWHAT=BWOSM*.5*(((A2+1)+(A2-1)*TAH)
     &        +(2*A1-A2-1)*EXP(-150*(LOMEGA-.25)**2))
C
C
C     **********  Lift Component  **********
C       Ref.   Y.Ikeda  et al  (JZK. No.143)
C
         IF( CM.LE.0.92 )THEN
           KAPA=0.0
         END IF
         IF( (CM.LE.0.97).AND.(CM.GT.0.92) )THEN
           KAPA=0.1
         END IF
         IF( CM.GT.0.97 )THEN
           KAPA=0.3
         END IF
         KN=6.28319*D/L+KAPA*(4.1*B/L-0.045)
         OG=OGD*D
         LO=0.3*D
         LR=0.5*D
         BLHAT=L*D*KN*LO*LR*FUN(IFN)*0.5/(NABLA*B**2)
     &         *SQRT(0.5*L*B)*(1-1.4*OG/LR+0.7*OG**2/(LO*LR))
C
C
                IF( RDP.EQ.2)THEN
C
C      ******* Eddy Making Component for type barge ships *******
C       Ref.  Y.Ikeda  et al(Proc. Ship Motion, Wave Loads and
C                          Propaulsive Performance in a Seaway,1984)
C           Y.Ikeda, T.Fujiwara, T.Katayama (Proc ISOPE Vol.3, 1993)
C
         OG=OGD*D
         DO I=1,N
           BE2(I)=(2/3.141592)*DX(I)**4/NABLA/BX(I)**2*(B/2/9.8)**.5
     &            *(HO(I)**2+1-OG/DX(I))*(HO(I)**2
     &            +(1-OG/DX(I))**2)*THTA*OMEGA
         END DO
         BEHAT=0
         DO I=1,N-1
           XHA=(X(I+1)-X(I))
           SBE=XHA*(BE2(I)+BE2(I+1))/2
           BEHAT=BEHAT+SBE
         END DO
         BEHAT=BEHAT*L/(X(N)-X(1))
         IF (FUN(IFN).EQ.0) then
           BEHAT=BEHAT
           GOTO 5215
         END IF
         AK=OMEGA/FUN(IFN)*SQRT(L/9.8)
         BEHAT=BEHAT*(.04*AK)**2/((.04*AK)**2+1)
 5215    CONTINUE
C
C
                  GOTO 5064
                END IF
C
C     *********** Eddy Making Component **********
C     Ref.   Y.Ikeda  et al  (JZK. No.142 , JZK. No.143)
C     Ref.   T.Katayama et al  (PJASNAOE No.33)
C
         DO J=1 , N
           AHO=HO(J)/(1-OGD)
           SGUMA=(SIG(J)-OGD)/(1-OGD)
           E=(AHO-1)/(AHO+1)
           E2=E**2
           A=4*SGUMA*(1-E2)/3.1415+E2
           O=-A/(A+3)
           O2=SQRT(O**2-(A-1)/(A+3))
           A3=O+O2
           A1=E*(1+A3)
           AM=BX(J)/(1+A1+A3)*0.5
           AA1=A1*(1+A3)/A3*0.25
           IF( AA1.GT.1 )THEN
             AA1=1
           END IF
           IF( AA1.LE.-1 )THEN
             AA1=-1
           END IF
           DO I=1 , 2
             CALL ARC(AA1,ACS)
             LTHETA=.5*ACS
             IF( I.EQ.1 )THEN
               LTHETA=0
             END IF
             AH=1+A1**2+9*A3**2+2*A1*(1-3*A3)*COS(2*LTHETA)
     &          -6*A3*COS(4*LTHETA)
             AA=-2*A3*COS(5*LTHETA)+A1*(1-A3)*COS(3*LTHETA)
     &          +((6-3*A1)*A3**2+(A1**2-3*A1)*A3+A1**2)*COS(LTHETA)
             BB=-2*A3*SIN(5*LTHETA)+A1*(1-A3)*SIN(3*LTHETA)
     &          +((6+3*A1)*A3**2+(3*A1+A1**2)*A3+A1**2)*SIN(LTHETA)
             V(I)=2*AM*SQRT(AA**2+BB**2)/AH
             RMAX1(I)=AM*SQRT(((1+A1)*SIN(LTHETA)-A3*SIN(3*LTHETA))**2
     &                +((1-A1)*COS(LTHETA)+A3*COS(3*LTHETA))**2)
           END DO
           RMAX=RMAX1(1)
           VMAX=V(1)
           IF( RMAX1(1).LE.RMAX1(2) )THEN
             GOTO 5159
           END IF
           GOTO 5161
 5159      RMAX=RMAX1(2)
           VMAX=V(2)
 5161      RMEAN=2*DX(J)*(1-OGD)*SQRT(AHO*SGUMA/3.1415)
           P1=VMAX/RMEAN
           P2=RMAX/RMEAN
           PP3=P1+P2
           GAMMA=(1+4*EXP(-165000*(1-SGUMA)**2))*PP3
           CP=0.5*(0.87*EXP(-GAMMA)-4*EXP(-0.187*GAMMA)+3)
           SIT=20*(SIG(J)-0.7)
           CALL TANH(TAH,SIT)
           F1=0.5*(1+TAH)
           F2=0.5*(1-COS(3.1415*SIG(J)))-1.5*(1-EXP(-5*(1-SIG(J))))
     &        *SIN(3.1415*SIG(J))**2
           IF( SIG(J).GT.1 )THEN
             SIG(J)=1.
           END IF
           R=2*DX(J)*SQRT(HO(J)*(SIG(J)-1)/(-0.8584))
           RD=R/DX(J)
           IF( (HO(J).LE.1).AND.(RD.GE.AHO) )THEN
             R=.5*BX(J)
           END IF
           IF( (HO(J).GE.1).AND.(RD.GE.1) )THEN
             R=DX(J)
           END IF
           RD=R/DX(J)
           CR1(J)=RMAX**2/DX(J)**2*CP*((1-F1*RD)*(1-3/2*OGD-F1*RD)
     &           +F2*(HO(J)-F1*RD)**2)
           IF( CR1(J).LT.0 )THEN
             CR1(J)=0
           END IF
         END DO
         DO K=1 , 21
           X1(K)=0+.5*(K-1)
           MAX=N
           MAXX=SC
           NN=MAX
           M1=1
           M2=0
           XXX=X1(K)
           DO JJ=1 , MAXX
             X11(JJ)=X(JJ)
             Y11(JJ)=CR1(JJ)
           END DO
           CALL HOKAN1(NN,IL,XXX,YX,Y2,M1,M2,X11,Y11,WX,WY)
           CR(K)=Y2
           DAM=YX
         END DO
         CR(1)=1.5*(1-3/2*OGD)
         CR(21)=1.5*(1-3/2*OGD)
         SAM=0
         DO K=1,10
           K2=2*K
           K1=K2-1
           K3=K2+1
           SAM1=CR(K1)+4*CR(K2)+CR(K3)
           SAM=SAM+SAM1
         END DO
         CRT=SAM/60
         BEHAT=4*L*D**4/3/3.1415*OMEGA*SQRT(B/19.6)/NABLA/B**2*CRT*THTA
         DO ISSN=1,N
           BE2(ISSN)=CR1(ISSN)*4*D**4/3/3.1415*OMEGA
     &               *SQRT(B/19.6)/NABLA/B**2*THTA
         END DO
C
         IF (FUN(IFN).EQ.0) then
           BEHAT=BEHAT
           GOTO 5211
         END IF
C
         AK=OMEGA/FUN(IFN)*SQRT(L/9.8)
         BEHAT=BEHAT*(.04*AK)**2/((.04*AK)**2+1)
         DO ISSN=1,N
           BE2(ISSN)=BE2(ISSN)*(.04*AK)**2/((.04*AK)**2+1)
         END DO
 5211    CONTINUE
C
C
 5064           IF( BBK.LE.1E-10)THEN
                  GOTO 5067
                END IF
                IF( TYPE_BK.EQ.1)THEN
C
C      ********** Damping due to Bilge Keels **********
C      Ref.   Y.Ikeda  et al  (JZK. No.161 , JZK. No.165)
C
         XBK(1)=XBK1
         XBK(11)=XBK2
         DO I=2,10
           XBK(I)=XBK(I-1)+(XBK2-XBK1)*.1
         END DO
         MAX=N
         DO I=1,11
           MAXX=SC
           NN=MAX
           M1=1
           M2=0
           XXX=XBK(I)

           DO JJ=1,MAXX
             X11(JJ)=X(JJ)
             Y11(JJ)=HO(JJ)
           END DO
           CALL HOKAN1(NN,IL,XXX,YX,Y2,M1,M2,X11,Y11,WX,WY)
           HO1=Y2
           DAM=YX
           MAXX=SC
           NN=MAX
           M1=1
           M2=0
           XXX=XBK(I)
           DO JJ=1,MAXX
             X11(JJ)=X(JJ)
             Y11(JJ)=SIG(JJ)
           END DO
           CALL HOKAN1(NN,IL,XXX,YX,Y2,M1,M2,X11,Y11,WX,WY)
           SG1=Y2
           DAM=YX
           MAXX=SC
           NN=MAX
           M1=1
           M2=0
           XXX=XBK(I)
           DO JJ=1,MAXX
             X11(JJ)=X(JJ)
             Y11(JJ)=DX(JJ)
           END DO
           CALL HOKAN1(NN,IL,XXX,YX,Y2,M1,M2,X11,Y11,WX,WY)
           DX1=Y2
           DAM=YX
           MAXX=SC
           NN=MAX
           XXX=XBK(I)
           M1=1
           M2=0
           XXX=XBK(I)
           DO JJ=1,MAXX
             X11(JJ)=X(JJ)
             Y11(JJ)=BX(JJ)
           END DO
           CALL HOKAN1(NN,IL,XXX,YX,Y2,M1,M2,X11,Y11,WX,WY)
           BX1=Y2
           DAM=YX

           R=2*DX1*sqrt(HO1*(SG1-1)/(-.8585))
           RD=R/DX1
           IF((HO1.LE.1).AND.(RD.GE.HO1)) then
             R=.5*BX1
           END IF
           IF((HO1.GT.1).AND.(RD.GE.1)) then
             R=DX1
           END IF
           RD=R/DX1
           F=1+.3*EXP(-160*(1-SG1))
           RBK=DX1*SQRT((HO1-.2929*RD)**2+(1-OGD-.2929*RD)**2)
           M1=RD
           M2=OGD
           M3=1-M1-M2
           M4=HO1-M1
           M5=(.414*HO1+.0651*M1**2-(.382*HO1+.0106)*M1)
     &        /((HO1-.215*M1)*(1-.215*M1))
           M6=(.414*HO1+.0651*M1**2-(.382+.0106*HO1)*M1)
     &        /((HO1-.215*M1)*(1-.215*M1))
           SO=.3*(3.1415*F*RBK*THTA)+1.95*BBK
           M7=SO/DX1-.25*3.1415*M1
           R1=.25*3.1415*R
           IF(SO.LT.R1) then
              M7=0
           END IF
           M8=M7+.414*M1
           IF(SO.LT.R1) then
              M8=M7+1.414*(1-cos(SO/R))*M1
           END IF
           A=(M3+M4)*M8-M7**2
           BB=M4**3/3/(HO1-.215*M1)+(1-M1)**2*(2*M3-M2)
     &        /6/(1-.215*M1)+M1*(M3*M5+M4*M6)
           CPPLAS=1.2
           CPMINS=-22.5*BBK/(3.1415*RBK*F*THTA)-1.2
           CD=CPPLAS-CPMINS
C
C      *** BBKhat for unit length ***
C
           RATIO(I)=RBK*BBK*CD/(RBK*BBK*CD+.5*DX1**2
     &            *(-A*CPMINS+BB*CPPLAS))
           BBKHAT(I)=8*RBK**2*OMEGA*SQRT(B/19.6)
     &               *THTA*F**2/(3*3.1415*NABLA*B**2)
     &               *(RBK*BBK*CD+.5*DX1**2*(-A*CPMINS+BB*CPPLAS))
         END DO
C
C      *** BBKhat for Three Dimentional Ship Form ***
C
         SAM=0
         DO I=1,5
           I2=2*I
           SAM1=BBKHAT(I2-1)+4*BBKHAT(I2)+BBKHAT(I2+1)
           SAM=SAM+SAM1
         END DO
         BKHAT=SAM*(XBK2-XBK1)*.1/3*L*.1
C
C
                ELSE
C
C          ********** Damping due to Bilge Keels **********
C          Ref. Y.Ikeda, T.Katayama et al  (JKSNAJ. No.222)
C
         BBK0=0
         BBK1=0
         BBK2=0
         BBK3=0
         RO=102
         KG=(1-OGD)*D
         DO P=SSBK1,SSBK2
           SGUMA=(SIG(P)-OGD)/(1-OGD)
C      ******* normal pressure component *******
C      ******* F1,R,CD----->FBN
           K=NUBK(P)
           XX0=XX(P,K)
           XX1=BKB(P)
           XX2=XX(P,K+1)
           YY0=YY(P,K)
           YY1=BKT(P)
           YY2=YY(P,K+1)
             IF(XX2.EQ.XX1) THEN
               XX2=XX(P,K+2)
               YY2=YY(P,K+2)
             END IF
             LOW=102
             F1=1+.3*EXP((-1)*160*(1-SGUMA))
             R=SQRT((KG-YY1)**2+XX1**2)
             CD=22.5*BBK/(3.14159*F1*R*THTA)+2.4
             FBN=8/3*RO*R**2*BBK*THTA*F1**2*CD*OMEGA/6.28
C       ******* XA,YA ***** Quadric approximation of a hull surface inclination
             IF(XX0.EQ.XX2) THEN
               XA=0
               YA=YY1
               GOTO 3550
             END IF
             IF(YY0.EQ.YY2) THEN
               XA=XX1
               YA=KG
               GOTO 3550
             END IF
             FD0=YY0/((XX0-XX1)*(XX0-XX2))
             FD1=YY1/((XX1-XX0)*(XX1-XX2))
             FD2=YY2/((XX2-XX0)*(XX2-XX1))
             FD=2*XX1*(FD0+FD1+FD2)-(XX1+XX2)
     &          *FD0-(XX0+XX2)*FD1-(XX0+XX1)*FD2
             XA=(XX1+FD*(YY1-KG))/(1+FD**2)
             YA=(FD*XX1+FD**2*(YY1-KG)+KG)/(1+FD**2)
C        ******* LB1,BN
 3550        LB1=SQRT((XA-XX1)**2+(YA-YY1)**2)
             BN=FBN*LB1*2
C
C        ******* negative pressure component *******
C        ******* S0,CPN----->FBSN
             S0=.3*3.14159*F1*R*THTA+1.95*BBK
             CPN=1.2-CD
             FBSN=8/3*RO*R**2*F1**2*THTA*CPN*OMEGA/6.28
C        ******* starboard side
C        ******* XX0,YY0,XX1,YY1
             MPN=0
             LB2S0=0
             OWARI=0
             BSN=0
             BSNS=0
             BSNP=0
             K=NUBK(P)
             XX0=BKB(P)
             YY0=BKT(P)
             XX1=XX(P,K+1)
             YY1=YY(P,K+1)
             GOTO 3950
 3900        K=K+1
             XX0=XX(P,K)
             YY0=YY(P,K)
             XX1=XX(P,K+1)
             YY1=YY(P,K+1)
             GOTO 3950
 3920        XX1=(XX1-XX0)*(S0+LB2-LB2S0)/LB2+XX0
             YY1=(YY1-YY0)*(S0+LB2-LB2S0)/LB2+YY0
C         ******* XM,YM,XA,YA
 3950        XM=(XX0+XX1)/2
             YM=(YY0+YY1)/2
             IF(XX0.EQ.XX1) THEN
               XA=0
               YA=YM
               GOTO 4010
             END IF
             IF(YY0.EQ.YY1) THEN
               XA=XM
               YA=KG
               GOTO 4010
             END IF
             XA=(YM+(XX0-XX1)/(YY0-YY1)*XM-KG)*(XX0-XX1)
     &          *(YY0-YY1)/((YY0-YY1)**2+(XX0-XX1)**2)
             YA=(YY0-YY1)/(XX0-XX1)*XA+KG
C         ******* LB2,LB3----->BSNS
 4010        LB2=SQRT((XX0-XX1)**2+(YY0-YY1)**2)
             LB2S0=LB2S0+LB2
             IF((OWARI.EQ.0).AND.(S0.LT.LB2S0)) THEN
               OWARI=1
               GOTO 3920
             END IF
             LB3=SQRT((XM-XA)**2+(YM-YA)**2)
             BSNS=BSNS+FBSN*LB2/LB3*(XA*(YM-YA)-(YA-KG)*(XM-XA))
             IF(OWARI.EQ.1) THEN
               OWARI=0
               GOTO 4110
             END IF
             GOTO 3900
C         ******* port side
C         ******* XX0,YY0,XX1,YY1
 4110        LB2S0=0
             K=2*NUMPT(P)-1-NUBK(P)
             XX0=-BKB(P)
             YY0=BKT(P)
             XX1=XX(P,K+1)
             YY1=YY(P,K+1)
             GOTO 4160
 4140        K=K+1
             XX0=XX(P,K)
             YY0=YY(P,K)
             XX1=XX(P,K+1)
             YY1=YY(P,K+1)
 4160        IF(YY1.GT.D) THEN
               OWARI=1
               XX1=(XX1-XX0)*(D-YY0)/(YY1-YY0)+XX0
               YY1=D
             END IF
             GOTO 4210
 4180        XX1=(XX1-XX0)*(S0+LB2-LB2S0)/LB2+XX0
             YY1=(YY1-YY0)*(S0+LB2-LB2S0)/LB2+YY0
C         ******* XM,YM,XA,YA
 4210         XM=(XX0+XX1)/2
             YM=(YY0+YY1)/2
             IF(XX0.EQ.XX1) THEN
               XA=0
               YA=YM
               GOTO 4270
             END IF
             IF(YY0.EQ.YY1) THEN
               XA=XM
               YA=KG
               GOTO 4270
             END IF
             XA=(YM+(XX0-XX1)/(YY0-YY1)*XM-KG)*(XX0-XX1)
     &          *(YY0-YY1)/((YY0-YY1)**2+(XX0-XX1)**2)
             YA=(YY0-YY1)/(XX0-XX1)*XA+KG
C         ******* LB2,LB3----->BSNP
 4270        LB2=SQRT((XX0-XX1)**2+(YY0-YY1)**2)
             LB2S0=LB2S0+LB2
             IF((OWARI.EQ.0).AND.(S0.LT.LB2S0)) THEN
               OWARI=1
               GOTO 4180
             END IF
             LB3=SQRT((XM-XA)**2+(YM-YA)**2)
             BSNP=BSNP+FBSN*LB2/LB3*(XA*(YM-YA)-(YA-KG)*(XM-XA))
             IF(OWARI.EQ.1) THEN
               OWARI=0
               GOTO 4350
             END IF
           GOTO 4140
 4350      BSN=BSNS+BSNP
C
C      ******* positive pressure component *******
C      ******* CPP
           CPP=1.2
C      ******* starboard side
C      ******* GLWBK
           GLWBK=0
           GLBKK0=0
           GLBKK1=0
           BSP=0
           BSPS=0
           BSPP=0
           DO K=1,NUBK(P)
             XX0=XX(P,K)
             YY0=YY(P,K)
             XX1=XX(P,K+1)
             YY1=YY(P,K+1)
             IF(YY1.GE.D) THEN
               goto 4520
             END IF
             IF((YY1.LT.D).AND.(YY0.GE.D)) then
               XX0=(XX1-XX0)*(D-YY0)/(YY1-YY0)+XX0
               YY0=D
             END IF
C
             IF(K.EQ.NUBK(P)) THEN
               XX1=BKB(P)
               YY1=BKT(P)
             END IF
             GLWBK=GLWBK+SQRT((XX0-XX1)**2+(YY0-YY1)**2)
 4520      END DO
C      ******* GLBKK0,GLBKK1,FBSP0,FBSP1
           DO K=1,NUBK(P)
             XX0=XX(P,K)
             YY0=YY(P,K)
             XX1=XX(P,K+1)
             YY1=YY(P,K+1)
             IF(YY1.GE.D) then
               goto 4740
             END IF
             IF ((YY1.LT.D).AND.(YY0.GE.D)) then
               XX0=(XX1-XX0)*(D-YY0)/(YY1-YY0)+XX0
               YY0=D
             END IF
             IF(K.EQ.NUBK(P)) THEN
               XX1=BKB(P)
               YY1=BKT(P)
             END IF
             GLBKK0=GLBKK1
             GLBKK1=GLBKK0+SQRT((XX0-XX1)**2+(YY0-YY1)**2)
             FBSP0=8/3*RO*R**2*THTA*F1**2*CPP*OMEGA/6.28*GLBKK0/GLWBK
             FBSP1=8/3*RO*R**2*THTA*F1**2*CPP*OMEGA/6.28*GLBKK1/GLWBK
C      ******* XM,YM,XA,YA,LB2,LB3----->BSP
             XM=(FBSP0+2*FBSP1)/(3*FBSP0+3*FBSP1)*(XX1-XX0)+XX0
             YM=(FBSP0+2*FBSP1)/(3*FBSP0+3*FBSP1)*(YY1-YY0)+YY0
             IF (XX0.EQ.XX1) THEN
               XA=0
               YA=YM
               GOTO 4700
             END IF
             IF(YY0.EQ.YY1) THEN
               XA=XM
               YA=KG
               GOTO 4700
             END IF
             XA=(YM+(XX0-XX1)/(YY0-YY1)*XM-KG)*(XX0-XX1)
     &          *(YY0-YY1)/((YY0-YY1)**2+(XX0-XX1)**2)
             YA=(YY0-YY1)/(XX0-XX1)*XA+KG
 4700        LB2=SQRT((XX0-XX1)**2+(YY0-YY1)**2)
             LB3=SQRT((XM-XA)**2+(YM-YA)**2)
             BSPS=BSPS+(FBSP0+FBSP1)/2*LB2/LB3
     &            *(XA*(YM-YA)-(YA-KG)*(XM-XA))
 4740      END DO
C       ******* port side
C       ******* GLWBK
           GLWBK=0
           GLBKK0=0
           GLBKK1=0
           OWARI=0


           KSTA=NUMPT(P)
           KEND=2*NUMPT(P)-1-NUBK(P)
           DO K=KSTA,KEND
             XX0=XX(P,K)
             YY0=YY(P,K)
             XX1=XX(P,K+1)
             YY1=YY(P,K+1)
             IF(K.EQ.KEND) THEN
               XX1=-BKB(P)
               YY1=BKT(P)
             END IF
             GLWBK=GLWBK+SQRT((XX0-XX1)**2+(YY0-YY1)**2)
           END DO
C       ******* GLBKK0,GLBKK1,FBSP0,FBSP1
           DO K=KSTA,KEND
             XX0=XX(P,K)
             YY0=YY(P,K)
             XX1=XX(P,K+1)
             YY1=YY(P,K+1)
             IF(K.EQ.KEND) THEN
               XX1=-BKB(P)
               YY1=BKT(P)
             END IF
             GLBKK0=GLBKK1
             GLBKK1=GLBKK0+SQRT((XX0-XX1)**2+(YY0-YY1)**2)
             FBSP0=8/3*RO*R**2*THTA*F1**2*CPP*OMEGA/6.28*GLBKK0/GLWBK
             FBSP1=8/3*RO*R**2*THTA*F1**2*CPP*OMEGA/6.28*GLBKK1/GLWBK
C       ******* XM,YM,XA,YA,LB2,LB3----->BSP
             XM=(FBSP0+2*FBSP1)/(3*FBSP0+3*FBSP1)*(XX1-XX0)+XX0
             YM=(FBSP0+2*FBSP1)/(3*FBSP0+3*FBSP1)*(YY1-YY0)+YY0
             IF (XX0.EQ.XX1) THEN
               XA=0
               YA=YM
               GOTO 4990
             END IF
             IF (YY0.EQ.YY1) THEN
               XA=XM
               YA=KG
               GOTO 4990
             END IF
             XA=(YM+(XX0-XX1)/(YY0-YY1)*XM-KG)*(XX0-XX1)
     &       *(YY0-YY1)/((YY0-YY1)**2+(XX0-XX1)**2)
             YA=(YY0-YY1)/(XX0-XX1)*XA+KG
 4990        LB2=SQRT((XX0-XX1)**2+(YY0-YY1)**2)
             LB3=SQRT((XM-XA)**2+(YM-YA)**2)
             BSPP=BSPP+(FBSP0+FBSP1)/2*LB2/LB3
     &             *(XA*(YM-YA)-(YA-KG)*(XM-XA))
           END DO
           BSP=BSPS+BSPP
C
C********************************************************************
           BBK0=BBK1
           BBK1=BN+BSN+BSP
           BBK2=BBK0+BBK1
           IF (P.EQ.SSBK1) then
             BBK3=BBK3+L*(CSD(P)-XBK1)/10*BBK1
             GOTO 5140
           END IF
           BBK3=BBK3+L*(CSD(P)-CSD(P-1))/10*BBK2/2
 5140      IF (P.EQ.SSBK2) THEN
             BBKDAMP=BBK3+L*(XBK2-CSD(P))/10*BBK1
           END IF
C       *******
         END DO
         BKHAT=BBKDAMP/(RO*NABLA*B**2)*SQRT(B/19.6)
C       *******
C=======================================
C
C
                END IF
                B2KHAT=BKHAT
 5067           IF( BBK.LE.1E-10)THEN
                  B2KHAT=0.
                END IF
                 BB44H=BFHAT+BWHAT+BLHAT+BEHAT+B2KHAT
                 B44 =BB44H*(RO*NABLA*B**2)/SQRT(B/19.6)
                 BF  =BFHAT *(RO*NABLA*B**2)/SQRT(B/19.6)
                 BW  =BWHAT *(RO*NABLA*B**2)/SQRT(B/19.6)
                 BL  =BLHAT *(RO*NABLA*B**2)/SQRT(B/19.6)
                 BE  =BEHAT *(RO*NABLA*B**2)/SQRT(B/19.6)
                 BBBK=B2KHAT*(RO*NABLA*B**2)/SQRT(B/19.6)
C                  WRITE(*,100)KAI(IKAI),BB44H,BFHAT
                  WRITE(2,100)THETA1(ITHETA),FUN(IFN),RAM(IRAM)
     &                  ,KAI(IKAI),BB44H,BFHAT,BWHAT,BLHAT,BEHAT
     &                  ,B2KHAT,B44,BF,BW,BL,BE,BBBK
  100             FORMAT( 16(F20.9,','))
c                ANSWER=0
c                DO ISSN=1 , N-1
c                  ANSWER=ANSWER+(X(ISSN+1)-X(ISSN))
c     &            *(BE2(ISSN)+BE2(ISSN+1))/2./10.*L
c                END DO
              END DO
            END DO
          END DO
        END DO
        CLOSE(2)

        STOP
        END
C
C
C
C       ************   Subroutine   *************
C
C      ********** Lagrange 3 Points Interpolation **********
      SUBROUTINE HOKAN1(NN,IL,XXX,YX,Y2,M1,M2,X11,Y11,WX,WY)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION NUE,L,KAI,L0,L1,L2,NABLA,KW,KAI1,LTHETA,MGA,KN
     &    ,LO,LR,KG,LB1,LB2,LB3,LB2S0,M1,M2,M3,M4,M5,M6,M7,M8
      INTEGER BWOSM,TYPE_BK,RDP,CSDMIN,SSBK1,SSBK2
      PARAMETER SC=200,PT=200,NNRA=70,NNFN=20,NNKAI=35,NNRD=30
C
      DIMENSION  X11(SC),Y11(SC),WX(SC),WY(SC)
C
C
         N1=NN-1
         DO IL=2,N1
           IF(XXX.LE.X11(IL)) THEN
             GOTO 5324
           END IF
         END DO
 5324    I1=IL-1
         IF (XXX.GT.X11(N1)) THEN
           I1=NN-2
         END IF
         I2=I1+2
         DO IL=I1,I2
           II=IL+1-I1
           WX(II)=X11(IL)
           WY(II)=Y11(IL)
         END DO
         IF (M1.NE.1) THEN
           GOTO 5334
         END IF
           CALL LAG3(XXX,Y2,X11,Y11,WX,WY)
 5334     IF (M2.NE.1) THEN
           GOTO 5336
         END IF
         YX=0
 5336  RETURN
       END
C
C
C
      SUBROUTINE LAG3(XXX,Y2,X11,Y11,WX,WY)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER SC=200,PT=200,NNRA=70,NNFN=20,NNKAI=35,NNRD=30
      DOUBLE PRECISION NUE,L,KAI,L0,L1,L2,NABLA,KW,KAI1,LTHETA,MGA,KN
     &    ,LO,LR,KG,LB1,LB2,LB3,LB2S0,M1,M2,M3,M4,M5,M6,M7,M8
      INTEGER BWOSM,TYPE_BK,RDP,CSDMIN,SSBK1,SSBK2
      DIMENSION  X11(SC),Y11(SC),WX(SC),WY(SC)
C
         Y2=0
         DO I3=1,3
           W=1
           Z=1
           DO J3=1,3
             IF(J3.EQ.I3) THEN
               GOTO 5350
             END IF
             W=W*(XXX-WX(J3))
             Z=Z*(WX(I3)-WX(J3))
 5350      END DO
           Y2=Y2+WY(I3)*W/Z
         END DO
       RETURN
       END
C
C
C      ********** arc cos **********
      SUBROUTINE ARC(AA1,ACS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
         IF(ABS(AA1).EQ.1) THEN
           ACS=0
           GOTO 5361
         END IF
         ACS=-ATAN(AA1/SQRT(-AA1*AA1+1))+3.14159/2
 5361  RETURN
       END
C
C
C      ********** hyperbolic tangent **********
      SUBROUTINE TANH(TAH,SIT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
         TAH=(EXP(SIT)-EXP(-SIT))/(EXP(SIT)+EXP(-SIT))
       RETURN
       END
C
C