!*********************************************
! Buckling pressure of steel pipe with stiffeners 
! under external ground water pressure
! (Technical Standards for gates and penstocks)
!*********************************************
module defpi
implicit none
real(8),parameter::pi=3.14159265358979323846D0
end module defpi

program f90_stiffener
    use defpi
    implicit none
    real(8)::Es    !Elastic modulus of steel pipe
    real(8)::po    !Poisson's ratio
    real(8)::As    !Thermal expansion coefficient
    real(8)::DD    !Internal diameter of steel pipe
    real(8)::t0    !Design thickness of steel pipe
    real(8)::eps   !Margin thickness of steel pipe
    real(8)::Sigy  !Yield stress of steel pipe
    real(8)::Siga  !Allowable stress of steel pipe
    real(8)::Jc    !Joint efficiency
    real(8)::Temp  !Temperature change
    real(8)::Bg    !Plastic deformation modulus of rock
    real(8)::SigF  !Modified yield stress
    real(8)::t     !pipe thickness without margin
    real(8)::rm    !rm=(D0+t0)/2
    real(8)::r0d   !r0d=(D0+2*t0)/2
    real(8)::k0    !Gap size between outer surface of steel pipe and inner surface of backfilled concrete
    real(8)::Pk    !Buckling pressure of steel pipe
    real(8)::tr0   !Thickness of stiffner
    real(8)::tr    !Thickness of stiffner without corrosion allowance
    real(8)::hr    !Height of stiffner
    real(8)::ell   !Stiffener's interval
    integer::n     !Number of wrinkles
    real(8)::sig_cr !Backling stress of Stiffener
    real(8)::sig_cc !Average compressive stress of Stiffener
    real(8)::S0,beta,ais,TT,lambda,elld
    real(8)::coef,pd
    real(8)::cc,c1,c2
    character(len=50)::dummy

!HT100  Sigy=885.0D0  Siga=400.0D0
!SHY685 Sigy=685.0D0  Siga=330.0D0
!SM570  Sigy=450.0D0  Siga=240.0D0
!SM490  Sigy=315.0D0  Siga=175.0D0
!SM400  Sigy=235.0D0  Siga=130.0D0

    call getarg(1,dummy); read(dummy,*) DD
    call getarg(2,dummy); read(dummy,*) t0
    call getarg(3,dummy); read(dummy,*) Sigy
    call getarg(4,dummy); read(dummy,*) Siga
    call getarg(5,dummy); read(dummy,*) tr0
    call getarg(6,dummy); read(dummy,*) hr
    call getarg(7,dummy); read(dummy,*) ell

    Es=206000.0D0
    po=0.3D0
    As=1.2D-5
    eps=1.5D0
    Jc=1.00D0
    Temp=20.0D0
    Bg=1.0D0
    SigF=Sigy

    t=t0-eps
    rm=(DD+t0)/2.0D0
    r0d=(DD+2.0D0*t0)/2.0D0
    k0=(As*Temp+Bg*Siga*Jc/Es)*r0d/(1.0D0+Bg)
    !k0=0.4D-3*rm
    tr=tr0-eps
    sig_cr=SCR(k0,rm,Es,t,SigF,tr,hr)

    ! Critical water pressure
    beta=(3.0D0*(1.0D0-po*po))**0.25D0/sqrt(rm*t)
    S0=tr*(t+hr)
    ais=tr/12.0D0*(hr+t)**3.0D0
    cc=tr+1.56D0*sqrt(rm*t)
    c1=r0d*r0d/t-cc*r0d*r0d/(S0+1.56D0*t*sqrt(rm*t))
    c2=3.0D0/(3.0D0*(1.0D0-po*po))**0.75D0*(r0d/t)*1.5D0
    c2=c2*(sinh(beta*ell)+sin(beta*ell))/(cosh(beta*ell)-cos(beta*ell))
    c2=c2+2.0D0*r0d*r0d/(S0+1.56D0*t*sqrt(rm*t))
    coef=1.0D0/cc*(cc+2.0D0*c1/c2)
    TT=2.0D0/cc*c1/c2
    lambda=1.0D0-(1.0D0+TT)*(1.0D0+tr/1.56D0/sqrt(rm*t))/(1.0D0+S0/t/1.56D0/sqrt(rm*t))
    elld=(ell+1.56D0*sqrt(rm*t)*acos(lambda))*(1.0D0+0.037D0*sqrt(rm*t)/(ell+1.56D0*sqrt(rm*t)*acos(lambda))*t**3.0D0/ais)
    cc=(1.0D0-po*po)*r0d/Es/t

    write(6,'(4f10.3)') DD,t0,Sigy,Siga
    write(6,'(5f10.3)') tr0,hr,ell,elld,sig_cr

    do n=2,20
        c1=(1.0D0-po*po)/dble(n*n-1)/(1.0D0+dble(n*n)*elld*elld/pi/pi/r0d/r0d)**2.0D0
        c2=t*t/12.0D0/r0d/r0d*(dble(n*n-1)+(dble(2*n*n-1)-po)/(1.0D0+dble(n*n)*elld*elld/pi/pi/r0d/r0d))
        Pk=(c1+c2)/cc
        pd=Pk*coef/1.5D0
        sig_cc=pd*r0d*(tr+1.56D0*sqrt(rm*t))/(S0+1.56D0*t*sqrt(rm*t))
        write(6,'(i5,4f10.3)') n,Pk,pd,sig_cc,sig_cr/sig_cc
    end do

stop

contains

    real(8) function SCR(k0,rm,Es,t,SigF,tr,hr)
        real(8),intent(in)::k0,rm,Es,t,SigF,tr,hr
        real(8)::b,ai,ee,ii
        real(8)::f1,f2,fi,x1,x2,xi !Variables for bisection method
        real(8)::SigN  !Axial stress of stiffener at the moment of buckling

        b=2.0D0*0.78D0*sqrt(rm*t)+tr
        ee=(0.5D0*b*t*t+tr*hr*(t+0.5D0*hr))/(b*t+tr*hr)
        ai=b/3.0D0*((t-ee)**3.0D0-ee**3.0D0)+tr/3.0D0*((hr+t-ee)**3.0D0-(t-ee)**3.0D0)
        ii=sqrt(ai/(b*t+tr*hr))
        !Bisection method
        x1=0.0D0
        x2=SigF
        do
            xi=0.5D0*(x1+x2)
            f1=FUNC(k0,rm,x1,Es,t,SigF,ee,ii)
            f2=FUNC(k0,rm,x2,Es,t,SigF,ee,ii)
            fi=FUNC(k0,rm,xi,Es,t,SigF,ee,ii)
            if(f1*fi<0.0)x2=xi
            if(fi*f2<0.0)x1=xi
            SigN=xi
            if(abs(x2-x1)<1.0D-6)exit
        end do
        SCR=SigN*(1.0D0-r0d/ee*(SigF-SigN)/(1.0D0+3.0D0/2.0D0*pi)/Es)
    end function SCR


    real(8) function FUNC(k0,rm,SigN,Es,t,SigF, ee,ii)
        real(8),intent(in)::k0,rm,SigN,Es,t,SigF,ee,ii
        real(8)::a,b
        a=(k0/rm+SigN/Es)*(1.0D0+rm*rm/ii/ii*SigN/Es)**1.5D0
        b=1.68D0*rm/ee*(SigF-SigN)/Es*(1.0D0-0.25D0*rm/ee*(SigF-SigN)/Es)
        FUNC=a-b
    end function FUNC

end program f90_stiffener
