program f90_RCmntD
    ! *****************************************************
    ! Stress of RC against Moment and tensile axial force
    ! *****************************************************
    implicit none
    real(8),parameter::nn=15.0D0
    real(8)::hh,bb,As_t,dd_t,As_c,dd_c,MO,AF
    real(8)::sig_cc,sig_sc,sig_st
    real(8)::x !Distance between Compressive edge and Neutral axis
    character(len=50)::dummy

    ! Data input
    call getarg(1,dummy);read(dummy,*) hh   !Height of section
    call getarg(2,dummy);read(dummy,*) bb   !Width of section
    call getarg(3,dummy);read(dummy,*) dd_t !Distance between Tensile edge and Tensile re-bar
    call getarg(4,dummy);read(dummy,*) dd_c !Distance between Compressive edge and Compressive re-bar
    call getarg(5,dummy);read(dummy,*) As_t !Section area of Tensile re-bar
    call getarg(6,dummy);read(dummy,*) As_c !Section area of Compressive re-bar
    call getarg(7,dummy);read(dummy,*) MO   !Moment (positive value)
    call getarg(8,dummy);read(dummy,*) AF   !Tensile axial force (positive value)

    sig_cc=0.0D0 !Compressive stress of Concrete
    sig_sc=0.0D0 !Compressive stress of Re-bar
    sig_st=0.0D0 !Tensile stress of Re-bar

    call SEC_DOUBLE(MO,AF,nn,hh,bb,As_t,dd_t,As_c,dd_c,x,sig_cc,sig_sc,sig_st)

    write(6,'(a)') 'hh,bb,dd_t,dd_c'
    write(6,'(e15.7,3(",",e15.7))') hh,bb,dd_t,dd_c
    write(6,'(a)') 'As_t,As_c,x'
    write(6,'(e15.7,2(",",e15.7))') As_t,As_c,x
    write(6,'(a)') 'MO,AF'
    write(6,'(e15.7,1(",",e15.7))') MO,AF
    write(6,'(a)') 'sig_cc,sig_sc,sig_st'
    write(6,'(e15.7,2(",",e15.7))') sig_cc,sig_sc,sig_st

stop

contains

    subroutine SEC_DOUBLE(MO,AF,nn,hh,bb,As_t,dd_t,As_c,dd_c,x,sig_cc,sig_sc,sig_st)
        real(8),intent(in)::MO,AF,nn,hh,bb,As_t,dd_t,As_c,dd_c
        real(8),intent(out)::x,sig_cc,sig_sc,sig_st
        real(8)::ee,cc,cd,c1,c2,c3,ff,df,work
        real(8),parameter::eps=1.0D-6

        ee=MO/AF
        cc=hh/2.0D0-dd_t
        cd=hh/2.0D0-dd_c
        c1=3.0D0*(hh/2.0D0+ee)
        c2=6.0D0*nn/bb*(As_t*(ee-cc)+As_c*(ee+cd))
        c3=6.0D0*nn/bb*(As_t*(cc+hh/2.0D0)*(ee-cc)+As_c*(hh/2.0D0-cd)*(ee+cd))

        x=0.5*cd
        do
            ff=x**3-c1*x**2-c2*x+c3
            df=3.0D0*x**2-c1*x-c2
            x=x-ff/df
            if(abs(ff/df)<eps)exit
        end do
        work=nn*As_t/x*(cc+hh/2.0D0-x)-nn*As_c/x*(cd-hh/2.0D0+x)-bb*x/2.0D0
        sig_cc=AF/work
        sig_sc=nn*sig_cc/x*(cd-hh/2.0D0+x)
        sig_st=nn*sig_cc/x*(cc+hh/2.0D0-x)
    end subroutine SEC_DOUBLE

end program f90_RCmntD