Module global
 implicit none
 real, parameter :: Total_Mass=200., p=1.2, ve = 800., A=1., Cd=0.8, g =9.81, alpha=0.5
 integer, parameter :: number_of_equations=3, number_of_steps=20000, itermax = 100
 real :: time, del_time, time_total=20.
 real, dimension(number_of_equations) :: Q
 real, dimension(number_of_equations),target :: T
 real, pointer :: height, velocity, mass
 real :: acc, mexit, mfuel, acceleration
end module global

Module ODE_Solver
use global
implicit none
contains
    subroutine RHS(time,T,rhs_value,Number_of_equations)
    implicit none
    integer::Number_of_equations
    real:: T(Number_of_equations),RHS_value(Number_of_equations)
    real:: time, C
    C=(p * Cd * A)/2.
       if (mass < (Total_Mass - mfuel))then
         mexit = 0.
       endif
    rhs_value(1) = -mexit
    rhs_value(2) = velocity
    rhs_value(3) = ((mexit * ve) - mass * g - ((C * abs(velocity)*velocity))) / mass
    acceleration = rhs_value(3)
end subroutine RHS

   subroutine RKG(del_time,time,T,Q,Number_of_equations)
        implicit none
        integer::Number_of_equations, I, J
        real:: T(Number_of_equations),Q(Number_of_equations),RHS_value(Number_of_equations),A(4),B(4),C(4),D(4)
        real:: del_time, time , R1, R2
        data A/0.5,0.2928932188134525, &
        &  1.707106781186547,     &
        &  0.1666666666666667/
        data B/2.0,1.0,1.0,2.0/
        data C/0.5,0.2928932188134525, &
        &  1.707106781186547,0.5/
        data D/0.0,0.5,0.,0.5/
            do J=1,4
            time =time+D(J)*del_time
            call RHS(time,T,rhs_value,Number_of_equations)
                do I=1,Number_of_equations
                  R1=del_time*rhs_value(I)
                  R2=A(J)*(R1-B(J)*Q(I))
                  T(I)=T(I)+R2
                  Q(I)=Q(I)+R2*3.-C(J)*R1
                enddo
            enddo
        return
   end subroutine RKG
end Module ODE_Solver

module secant_solver
     use ODE_Solver
     implicit none
     contains
       subroutine secant(icount,mfuel_1,mfuel_2,height_1,height_ref)
       implicit none
       real :: mfuel_1, mfuel_2, height_1, height_2, height_ref, del_mfuel, height_max1, height_max2
       integer :: icount, i
            mfuel = mfuel_1
            do i = 1, number_of_steps
                 call RKG(del_time,time,T,Q,number_of_equations)
                 if (velocity <= 0.) exit
                 height_max1 = height
                 write(1,*)time,height,velocity,acceleration,mass
            enddo
            height_1 = height_max1
            mfuel = mfuel_2
            height = 0.
            velocity = 0.
            mass = 200.
            time = 0.
            mexit = 3.
            do i = 1, number_of_steps
                 call RKG(del_time,time,T,Q,number_of_equations)
                 if (velocity <= 0.) exit
                 height_max2 = height
                 write(2,*)time,height,velocity,acceleration,mass
            enddo
            height_2 = height_max2
            do icount = 1, itermax
                 del_mfuel = (height_ref - height_1) * (mfuel_1 - mfuel_2) / (height_1 - height_2)
                 mfuel = mfuel + del_mfuel * alpha
                 print*, abs(height_1-height_ref), abs(height_2-height_ref), del_mfuel
                 if (abs(height_1 - height_ref) < abs(height_2 - height_ref)) then
                      mfuel_2 = mfuel_1
                      height_2 = height_1
                 endif
                 mfuel_1 = mfuel
                 height = 0.
                  velocity = 0.
                  mass = 200.
                  time = 0.
                  mexit = 3.
                  do i = 1, number_of_steps
                       call RKG(del_time,time,T,Q,number_of_equations)
                       if (velocity <= 0.) exit
                       height_max1 = height
                       write(3,*)time,height,velocity,acceleration,mass
                 enddo
                 height_1 = height_max1
                 if (abs(del_mfuel) < 2.e-3) exit
            enddo
       RETURN
       end subroutine secant
end module secant_solver

program rocket
use secant_solver
implicit none
integer :: iter
real ::  mfuel_1, mfuel_2,height_1,height_ref,C
     mass => T(1)
     height => T(2)
     velocity => T(3)
     time = 0.
     mexit = 3.
     del_time = time_total / real(number_of_steps)
     height_ref = 200.
     height = 0.
     velocity = 0.
     mass = 200.
     C = (p*Cd*A)/ 2.
     acceleration =((mexit * ve) - (mass * g) - (C * abs(velocity) * velocity)) / mass
     Q(:) = 0.
          open(1, file='rocket.dat')
          write(1,*)'Variables = "Time" "Height" "Velocity" "Acceleration" "Mass"'
          write(1,*)'Zone I=', number_of_steps+1, ' Datapacking = point'
          write(1,*)time,height,velocity,acceleration,mass
     mfuel_1 = 25.
     mfuel_2 = 35.
     call secant(iter,mfuel_1,mfuel_2,height_1,height_ref)
     mfuel = mfuel_1
print*, time,height_1,velocity,acceleration,mass,mfuel
stop
end program rocket