program f90_TODAPP
    implicit none
    character::strcom*50
    integer::ndata
    real(8),allocatable::xd(:),yd(:),out_toda(:),out_shenton(:)
    integer i
    character::fnameR*50,fnameW*50

    call getarg(1,fnameR)
    call getarg(2,fnameW)

    !Data input
    open(11,file=fnameR,status='old')
    read(11,'(a)') strcom
    read(11,*) ndata
    allocate(xd(1:ndata),yd(1:ndata))
    allocate(out_toda(1:ndata),out_shenton(1:ndata))
    do i=1,ndata
        read(11,*) xd(i),yd(i)
        yd(i)=yd(i)+0.5D0
    end do
    close(11)

    !calculation of probability
    do i=1,ndata
        out_shenton(i)=SHENTON(xd(i))
        out_toda(i)=TODA(yd(i))
    end do

    open(12, file=fnameW, status='replace')
    write(12,'(a)') trim(adjustl(strcom))
    write(12,'(a)') 'Input,Input,Output,Output'
    write(12,'(a)') '%-point,None-exceed prob.,Non-exceed prob(Shenton),%-point(Toda)'
    do i=1,ndata
        write(12,'(f9.6,",",f9.6,",",f9.6,",",f9.6)') xd(i),yd(i),out_shenton(i),out_toda(i)
    end do
    close(12)

stop

contains

    real(8) function TODA(rx)
        !% point for specified non-exceedance probability
        real(8),intent(in)::rx

        integer::i
        real(8)::ry,ay,uay,b(0:10),sum

        if(rx==0.5D0)then
            uay=0.0D0
        else
            if(rx<0.5D0)then
                ry=rx
            else
                ry=1.0D0-rx
            end if
            ay=-log(4.0D0*ry*(1.0D0-ry))
            b(0)= 1.570796288D0
            b(1)= 0.03706987906D0
            b(2)=-0.0008364353589D0
            b(3)=-0.0002250947176D0
            b(4)= 0.000006841218299D0
            b(5)= 0.000005824238515D0
            b(6)=-0.00000104527497D0
            b(7)= 0.00000008360937017D0
            b(8)=-0.000000003231081277D0
            b(9)= 0.00000000003657763036D0
            b(10)=0.0000000000006936233982D0
            sum=0.0D0
            do i=0,10
                sum=sum+b(i)*ay**dble(i)
            end do
            uay=sqrt(ay*sum)
            if(rx<0.5D0)uay=-uay
        end if
        TODA=uay
    end function TODA


    real(8) function SHENTON(u)
        !Non-exceedance probability for specified % point
        real(8),intent(in)::u
        integer::nt,ii
        real(8)::ue,af,qh,ag,ab,q

        nt=30
        if(u>=0)ue=u
        if(u<0)ue=-u
        af=ue*ue
        qh=0.0D0
        ag=exp(-0.5D0*af)*0.3989422804014327D0
        ab=-1.0D0
        do ii=nt,1,-1
            qh=dble(ii)*af/(2.0D0*dble(ii)+1.0D0+ab*qh)
            ab=-ab
        end do
        qh=0.5D0-ag*ue/(1.0D0-qh)
        if(u<0.0D0)qh=1.0D0-qh
        q=qh
        SHENTON=1.0D0-q
    end function SHENTON


end program f90_TODAPP