program f90_ELLIP
!--------------------------------------------------------------------
! subroutine CALCOEF()
!   This is a subroutine to obtain the elements of matrix
!   for the iterarive calculation. 
! subroutine STEGJ()
!   This is a subroutine of solver for simultaneous equations. 
!--------------------------------------------------------------------
    implicit none
    integer,parameter::itmax=100
    real(8),parameter::eps=1D-6
    real(8),parameter::pi=3.14159265358979323846D0
    integer::i
    real(8)::aa,bb,tt
    real(8)::coefA,coefB,coefC
    real(8)::ff,dfda,dfdb,dfdt
    real(8)::gg,dgda,dgdb,dgdt
    real(8)::hh,dhda,dhdb,dhdt
    real(8)::ary(1:3,1:4)
    real(8)::x1,y1,x2,y2,x3,y3,xm,ym
    character(len=50)::fnameR,dummy
    real(8)::LX,LY,dx,dy

    call getarg(1,fnameR)
    call getarg(2,dummy)
        read(dummy,*) LX
    call getarg(3,dummy)
        read(dummy,*) LY
    call getarg(4,dummy)
        read(dummy,*) dx
    call getarg(5,dummy)
        read(dummy,*) dy

    open(11,file=fnameR,status='old')
    read(11,*) xm,ym
    read(11,*) x1,y1
    read(11,*) x2,y2
    read(11,*) x3,y3
    close(11)

    x1=x1*LX/dx
    y1=y1*LY/dy
    x2=x2*LX/dx
    y2=y2*LY/dy
    x3=x3*LX/dx
    y3=y3*LY/dy

    ary(1,1)=x1*x1
    ary(1,2)=x1*y1
    ary(1,3)=y1*y1
    ary(2,1)=x2*x2
    ary(2,2)=x2*y2
    ary(2,3)=y2*y2
    ary(3,1)=x3*x3
    ary(3,2)=x3*y3
    ary(3,3)=y3*y3
    ary(1,4)=1.0D0
    ary(2,4)=1.0D0
    ary(3,4)=1.0D0
    call STEGJ(3,ary)
    coefA=ary(1,4)
    coefB=ary(2,4)
    coefC=ary(3,4)

    ! Initial value
    aa=sqrt(x1*x1+y1*y1)
    bb=sqrt(x2*x2+y2*y2)
    tt=atan(y1/x1)
    ! Iteration
    do i=1,itmax
        call CALCOEF(aa,bb,tt,coefA,coefB,coefC,ary,3)
        call STEGJ(3,ary)
        aa=aa+ary(1,4)
        bb=bb+ary(2,4)
        tt=tt+ary(3,4)
        if(abs(ary(1,4))<eps.and.abs(ary(2,4))<eps.and.abs(ary(2,4))<eps)exit
    end do

    write(6,'(5(e15.7))') xm,ym,tt*180.0D0/pi,2.0D0*aa,2.0D0*bb

stop

contains

    subroutine CALCOEF(aa,bb,tt,coefA,coefB,coefC,ary,nn)
        integer,intent(in)::nn
        real(8),intent(in)::aa,bb,tt,coefA,coefB,coefC
        real(8),intent(out)::ary(1:nn,1:nn+1)

        real(8)::ff,dfda,dfdb,dfdt,gg,dgda,dgdb,dgdt,hh,dhda,dhdb,dhdt

        ff=bb*bb*cos(tt)*cos(tt)+aa*aa*sin(tt)*sin(tt)-aa*aa*bb*bb*coefA
        gg=2.0D0*sin(tt)*cos(tt)*(bb*bb-aa*aa)-aa*aa*bb*bb*coefB
        hh=bb*bb*sin(tt)*sin(tt)+aa*aa*cos(tt)*cos(tt)-aa*aa*bb*bb*coefC

        dfda= 2.0D0*aa*(sin(tt)*sin(tt)-bb*bb*coefA)
        dfdb= 2.0D0*bb*(cos(tt)*cos(tt)-aa*aa*coefA)
        dfdt= 2.0D0*sin(tt)*cos(tt)*(aa*aa-bb*bb)

        dgda=-2.0D0*aa*(2.0D0*sin(tt)*cos(tt)+bb*bb*coefB)
        dgdb= 2.0D0*bb*(2.0D0*sin(tt)*cos(tt)-aa*aa*coefB)
        dgdt=-2.0D0*(cos(tt)*cos(tt)-sin(tt)*sin(tt))*(aa*aa-bb*bb)

        dhda= 2.0D0*aa*(cos(tt)*cos(tt)-bb*bb*coefC)
        dhdb= 2.0D0*bb*(sin(tt)*sin(tt)-aa*aa*coefC)
        dhdt=-2.0D0*sin(tt)*cos(tt)*(aa*aa-bb*bb)

        ary(1,1)=dfda
        ary(1,2)=dfdb
        ary(1,3)=dfdt
        ary(2,1)=dgda
        ary(2,2)=dgdb
        ary(2,3)=dgdt
        ary(3,1)=dhda
        ary(3,2)=dhdb
        ary(3,3)=dhdt
        ary(1,4)=-ff
        ary(2,4)=-gg
        ary(3,4)=-hh

    end subroutine CALCOEF


    subroutine STEGJ(n,ary)
        integer,intent(in)::n
        real(8),intent(inout)::ary(1:n,1:n+1)

        integer::i,j,k,s
        real(8)::p,d,max,dumy

        do k=1,n
            max=0.0D0
            s=k
            do j=k,n
                if(abs(ary(j,k))>max)then
                    max=abs(ary(j,k))
                    s=j
                end if
            end do
            do j=1,n+1
                dumy=ary(k,j)
                ary(k,j)=ary(s,j)
                ary(s,j)=dumy
            end do
            p=ary(k,k)
            do j=k,n+1
                ary(k,j)=ary(k,j)/p
            end do
            do i=1,n
                if(i/=k)then
                    d=ary(i,k)
                    do j=k,n+1
                        ary(i,j)=ary(i,j)-d*ary(k,j)
                    end do
                end if
            end do
        end do
    end subroutine STEGJ

end program f90_ELLIP