module defpi
    implicit none
    real(8),parameter::pi=3.14159265358979323846D0
end module defpi

program f90_SREG
!-----------------------------------------------------------
! subroutine SREGRESSION()
!     This is a subroutine for single regression analysis.
! subroutine COELLIPSE()
!     This is a subroutine to define the caracteristics of 
!     probability ellipse. 
!-----------------------------------------------------------
    use defpi
    implicit none
    integer,parameter::NDmax=1000
    integer::i,ndata,nt,io
    real(8)::datax(1:NDmax),datay(1:NDmax)
    real(8)::aa,bb,rr,xm,ym,sx,sy
    real(8)::pp,rmax,rmin,tmax,tmin,xmin,yxmn
    character(len=50)::fnameR
    character(len=20)::dummy
    character(len=5)::nlx,nly

    call getarg(1,fnameR)
    call getarg(2,dummy)
    read(dummy,*) pp

    open(11,file=fnameR,status='old')
        read(11,*) nlx,nly
        do i=1,NDmax
            read(11,*,iostat=io) datax(i),datay(i)
            if(io<0)exit
        end do
    close(11)
    ndata=i-1
    nt=ndata

    ! If 'nlx'includes 'L', x-values are converted into common logarithms.
    if(index(nlx,"L")/=0)then
        do i=1,ndata
            datax(i)=log10(datax(i))
        end do
    end if
    ! If 'nly'includes 'L', y-values are converted into common logarithms.
    if(index(nly,"L")/=0)then
        do i=1,ndata
            datay(i)=log10(datay(i))
        end do
    end if

    call SREGRESSION(ndata,datax,datay,aa,bb,rr,xm,ym,sx,sy,nt)
    call COELLIPSE(pp,rr,xm,ym,sx,sy,rmax,rmin,tmax,tmin,xmin,yxmn)

    write(6,'("y=aa*x+bb")')
    write(6,'("aa=",e15.7)') aa
    write(6,'("bb=",e15.7)') bb
    write(6,'("rr=",e15.7)') rr
    write(6,'("ndata=",i5)') ndata
    write(6,'("xm=",e15.7)') xm
    write(6,'("ym=",e15.7)') ym
    write(6,'("sx=",e15.7)') sx
    write(6,'("sy=",e15.7)') sy
    write(6,'("pp=",e15.7)') pp
    write(6,'("rmax=",e15.7)') rmax
    write(6,'("rmin=",e15.7)') rmin
    write(6,'("tmax=",e15.7)') tmax
    write(6,'("tmin=",e15.7)') tmin
    write(6,'("(xm_ym)",2(e15.7))') xm,ym
    write(6,'("(x1_y1)",2(e15.7))') rmax*cos(tmax/180.0D0*pi),rmax*sin(tmax/180.0D0*pi)
    write(6,'("(x2_y2)",2(e15.7))') rmin*cos(tmin/180.0D0*pi),rmin*sin(tmin/180.0D0*pi)
    write(6,'("(x3_y3)",2(e15.7))') xmin,yxmn

stop

contains

    subroutine SREGRESSION(ndata,datax,datay,aa,bb,rr,xm,ym,sx,sy,nt)
        integer,intent(in)::ndata,nt
        real(8),intent(in)::datax(1:nt),datay(1:nt)
        real(8),intent(out)::aa,bb,rr
        real(8),intent(out)::xm,ym,sx,sy

        integer::i
        real(8)::x1,y1,x2,y2,xy,c1,c2,c3
        real(8)::rad,theta,xx,yy
        character(len=50)::fnameW

        !回帰式：y=aa*x+bb
        x1=0.0D0
        y1=0.0D0
        x2=0.0D0
        xy=0.0D0
        do i=1,ndata
            x1=x1+datax(i)
            y1=y1+datay(i)
            x2=x2+datax(i)*datax(i)
            xy=xy+datax(i)*datay(i)
        end do
        xm=x1/dble(ndata)
        ym=y1/dble(ndata)
        aa=(dble(ndata)*xy-x1*y1)/(dble(ndata)*x2-x1*x1)
        bb=(x2*y1-x1*xy)/(dble(ndata)*x2-x1*x1)
        c1=0.0D0
        c2=0.0D0
        c3=0.0D0
        do i=1,ndata
            c1=c1+(datax(i)-xm)*(datay(i)-ym)
            c2=c2+(datax(i)-xm)*(datax(i)-xm)
            c3=c3+(datay(i)-ym)*(datay(i)-ym)
        end do
        rr=c1/sqrt(c2*c3)
        sx=sqrt(c2/dble(ndata-1))
        sy=sqrt(c3/dble(ndata-1))
    end subroutine SREGRESSION


    subroutine COELLIPSE(pp,rr,xm,ym,sx,sy,rmax,rmin,tmax,tmin,xmin,yxmn)
        real(8),intent(in)::pp,rr,xm,ym,sx,sy
        real(8),intent(out)::rmax,rmin,tmax,tmin,xmin,yxmn

        integer,parameter::nn=3600
        integer::i
        real(8)::rad,theta,xx,yy,phi,dd
        character(len=50)::fnameW

        rmax=0.0D0
        rmin=1.0D10
        tmax=0.0D0
        tmin=0.0D0
        xmin=1.0D10
        do i=0,nn
            theta=dble(i)/dble(nn)*2.0D0*pi
            rad=sqrt(-2.0D0*(1.0D0-rr*rr)*log(1.0D0-pp)/(1.0D0-2.0D0*rr*sin(theta)*cos(theta)))
            xx=xm+sx*rad*cos(theta)
            yy=ym+sy*rad*sin(theta)
            dd=sqrt((xx-xm)**2.0D0+(yy-ym)**2.0D0)
            phi=atan((yy-ym)/(xx-xm))/pi*180.0D0
            if(rmax<=dd)then
                rmax=dd
                tmax=phi
            end if
            if(dd<=rmin)then
                rmin=dd
                tmin=phi
            end if
            if(xx-xm<=xmin)then
                xmin=xx-xm
                yxmn=yy-ym
            end if
        end do

        fnameW='_coellipse.txt'
        open(12, file=fnameW, status='replace')
        do i=0,360
            theta=dble(i)/180.0D0*pi
            rad=sqrt(-2.0D0*(1.0D0-rr*rr)*log(1.0D0-pp)/(1.0D0-2.0D0*rr*sin(theta)*cos(theta)))
            xx=xm+sx*rad*cos(theta)
            yy=ym+sy*rad*sin(theta)
            write(12,*) xx,yy
        end do
        close(12)

    end subroutine COELLIPSE

end program f90_SREG