program f90_ams_DD
    !*********************************************
    ! Buckling pressure of steel pipe under external ground water pressure
    ! (Technical Standards for gates and penstocks)
    !*********************************************
    implicit none
    real(8)::Es    !Elastic modulus of steel pipe
    real(8)::pos   !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)::Ess   !Modified elastic modulus
    real(8)::ny    !Coefficient of support effect
    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)::SigN  !Axial stress of pipe at the moment of buckling
    real(8)::Pk    !Buckling pressure of steel pipe
    real(8)::Phi,Psi,Omega !Coefficients for Amstutz's formula
    real(8)::f1,f2,fi,t1,t2,ti,rm1,rm2,rmi !Variables for bisection method
    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,*) Pk
    call getarg(3,dummy); read(dummy,*) Sigy
    call getarg(4,dummy); read(dummy,*) Siga

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

    Ess=Es/(1.0D0-pos*pos)
    ny=1.5D0-0.5D0*1.0D0/(1.0D0+0.002D0*Es/Sigy)/(1.0D0+0.002D0*Es/Sigy)
    SigF=ny*Sigy/sqrt(1.0D0-pos+pos*pos)
    r0d=(DD+2.0D0*t0)/2.0D0
    k0=(As*Temp+Bg*Siga*Jc/Es)*r0d/(1.0D0+Bg)
    !k0=0.4D-3*rm
    !Amstutz
    Phi=1.68D0
    Psi=0.25D0
    Omega=0.175D0
    !Bisection method
    t1=10.0D0
    t2=100.0D0
    do
        ti=0.5D0*(t1+t2)
!        write(6,*) t1,t2,ti
        rm1=(DD+t1+eps)/2.0D0
        rm2=(DD+t2+eps)/2.0D0
        rmi=(DD+ti+eps)/2.0D0
        f1=FUNC1(Pk,t1,rm1,k0,Ess,SigF,Phi,Psi,Omega)
        f2=FUNC1(Pk,t2,rm2,k0,Ess,SigF,Phi,Psi,Omega)
        fi=FUNC1(Pk,ti,rmi,k0,Ess,SigF,Phi,Psi,Omega)
        if(f1*fi<0.0D0)t2=ti
        if(f2*fi<0.0D0)t1=ti
        if(abs(t2-t1)<1.0D-6)exit
    end do
    t0=dble(int(ti+eps+1.0))
    write(6,'(5f10.3)') DD,t0,Sigy,Siga,Pk

stop

contains

    real(8) function FUNC1(Pk,t,rm,k0,Ess,SigF,Phi,Psi,Omega)
        real(8),intent(in)::Pk,t,rm,k0,Ess,SigF,Phi,Psi,Omega
        real(8)::x1,x2,xi,f1,f2,fi,SigN
        x2=SigF
        x1=x2
        do
            x1=x1-5.0
            f1=FUNC2(k0,rm,x1,Ess,t,SigF,Phi,Psi)
            f2=FUNC2(k0,rm,x2,Ess,t,SigF,Phi,Psi)
            if(f1*f2<0.0D0)exit
        end do
        if(x1<0.0D0)then
            write(6,*) 'Something is wrong!!'
            stop
        end if
        do
            xi=0.5D0*(x1+x2)
            f1=FUNC2(k0,rm,x1,Ess,t,SigF,Phi,Psi)
            f2=FUNC2(k0,rm,x2,Ess,t,SigF,Phi,Psi)
            fi=FUNC2(k0,rm,xi,Ess,t,SigF,Phi,Psi)
            if(f1*fi<0.0D0)x2=xi
            if(fi*f2<0.0D0)x1=xi
            if(abs(x2-x1)<1.0D-6)exit
        end do
        SigN=xi
        FUNC1=Pk-SigN/(rm/t)/(1.0D0+Omega*2.0D0*rm/t*(SigF-SigN)/Ess)
    end function FUNC1

    real(8) function FUNC2(k0,rm,SigN,Ess,t,SigF,Phi,Psi)
        real(8),intent(in)::k0,rm,SigN,Ess,t,SigF,Phi,Psi
        real(8)::a,b
        a=(k0/rm+SigN/Ess)*(1.0D0+12.0D0*rm*rm/t/t*SigN/Ess)**1.5D0
        b=Phi*2.0D0*rm/t*(SigF-SigN)/Ess*(1.0D0-Psi*2.0D0*rm/t*(SigF-SigN)/Ess)
        FUNC2=a-b
    end function FUNC2

end program f90_ams_DD
