!------------------!
 module systemparam
 implicit none

 real(8), parameter :: pi=3.141592653589793d0 
 real(8), parameter :: gm=3.987d14  ! Earth's mass times constant of gravity
 real(8), parameter :: arocket=5.d0 ! deceleration due to rocket motor
 real(8), parameter :: dragc=8.d-4  ! Air drag coefficient
 real(8), parameter :: re=6.378d6   ! Earth's radius

 real(8) :: dt,dt2  ! time step, half of the time step
 real(8) :: tbrake  ! run-time of rocket motor

 end module systemparam
!----------------------!


!======================!
 program satellitecrash
!------------------------------------------------------------------!
!calculates the trajectory of a satellite that is initially in a   !
!circular orbit and begins to brake at the start of the calculation!
!------------------------------------------------------------------!
 use systemparam
 implicit none

 integer :: i,nstp,wstp
 real(8) :: a,r,t,x,y,vx,vy,r0,tmax

 print*,'Initial altitude of satellite (km)'; read*,r0
 r0=r0*1.d3+re
 print*,'Rocket motor run-time (seconds)'; read*,tbrake
 print*,'Time step delta-t for RK integration (seconds)';read*,dt
 dt2=dt/2.d0
 print*,'Writing results every Nth step; give N';read*,wstp
 print*,'Maximum integration time (hours)';read*,tmax
 tmax=tmax*3600.d0

 x=r0
 y=0.d0
 vx=0.d0
 vy=sqrt(gm/r0)
 nstp=int(tmax/dt)

 open(1,file='sat.dat',status='replace')
 do i=0,nstp
   call polarposition(x,y,r,a)
   if (r > re) then
     t=dble(i)*dt
     if (mod(i,wstp)==0) write(1,1)t,a,(r-re)/1.d3,sqrt(vx**2+vy**2)
     1 format(f12.3,'  ',f12.8,'  ',f14.6,'  ',f12.4)
     call rkstep(t,x,y,vx,vy)
   else
     print*,'The satellite has successfully crashed!'
     goto 2
   end if
 end do
 print*,'The satellite did not crash within the specified time.'
 2 close(1)

 end program satellitecrash
!==========================!

!-----------------------------------!
 subroutine rkstep(t0,x0,y0,vx0,vy0)
!---------------------------------------!
!Integrates the equations of motion one !
!time step, using the Runge-Kutta method!
!---------------------------------------!
 use systemparam
 implicit none

 real(8) :: x0,y0,x1,y1,vx0,vy0,vx1,vy1,t0,th,t1,ax,ay
 real(8) :: kx1,kx2,kx3,kx4,ky1,ky2,ky3,ky4
 real(8) :: lx1,lx2,lx3,lx4,ly1,ly2,ly3,ly4

 t1=t0+dt; th=t0+dt2   
 call accel(x0,y0,vx0,vy0,t0,ax,ay) 
 kx1=dt2*ax 
 ky1=dt2*ay
 lx1=dt2*vx0 
 ly1=dt2*vy0
 call accel(x0+lx1,y0+ly1,vx0+kx1,vy0+ky1,th,ax,ay)
 kx2=dt2*ax; ky2=dt2*ay
 lx2=dt2*(vx0+kx1)
 ly2=dt2*(vy0+ky1)
 call accel(x0+lx2,y0+ly2,vx0+kx2,vy0+ky2,th,ax,ay)
 kx3=dt*ax
 ky3=dt*ay
 lx3=dt*(vx0+kx2)
 ly3=dt*(vy0+ky2)
 call accel(x0+lx3,y0+ly3,vx0+kx3,vy0+ky3,t1,ax,ay)
 kx4=dt2*ax
 ky4=dt2*ay
 lx4=dt2*(vx0+kx3)
 ly4=dt2*(vy0+ky3)
 x1=x0+(lx1+2.d0*lx2+lx3+lx4)/3.d0
 y1=y0+(ly1+2.d0*ly2+ly3+ly4)/3.d0
 vx1=vx0+(kx1+2.d0*kx2+kx3+kx4)/3.d0
 vy1=vy0+(ky1+2.d0*ky2+ky3+ky4)/3.d0

 x0=x1; y0=y1; vx0=vx1; vy0=vy1

 end subroutine rkstep
!---------------------!

!-----------------------------------!
 subroutine accel(x,y,vx,vy,t,ax,ay)
!------------------------------------------------------!
!calculates the x and y components of the acceleration,!
!due to gravitation, air drag, and brake-motor thrust  !
!------------------------------------------------------!
 use systemparam
 implicit none

 real(8) :: x,y,vx,vy,t,ax,ay,r,v1,v2,r3,ad,airdens

 r=sqrt(x**2+y**2)       
 v2=vx**2+vy**2
 v1=sqrt(v2)

   !*** evaluates the acceleration due to gravitation
 r3=1.d0/r**3
 ax=-gm*x*r3           
 ay=-gm*y*r3  

   !*** evaluates the acceleration due to air drag
 if (v1 > 1.d-12) then
   ad=dragc*airdens(r)*v2
   ax=ax-ad*vx/v1
   ay=ay-ad*vy/v1
 end if

   !*** evaluates the acceleration due to rocket motor thrust
 if (t < tbrake .and. v1 > 1.d-12) then
   ax=ax-arocket*vx/v1
   ay=ay-arocket*vy/v1
 end if

 end subroutine accel
!--------------------!

!---------------------------!
 real(8) function airdens(r)
!-------------------------------------------------!
!the assumed air-density as a function of altitude!
!-------------------------------------------------!
 use systemparam
 implicit none

 real(8), parameter :: k1=1.2d4
 real(8), parameter :: k2=2.2d4

 real(8) :: r,rr

 if (r < re) return
 airdens=1.225d0*exp(-((r-re)/k1+((r-re)/k2)**1.5))

 end function airdens
!--------------------!

!---------------------------------!
 subroutine polarposition(x,y,r,a)
!------------------------------------------------------!
!converts the x and y coordinates to a distance r from !
!the center of the earth and a normalized angle a (0-1)!
!------------------------------------------------------!
 use systemparam
 implicit none
 
 real(8) :: x,y,r,a

 r=sqrt(x**2+y**2)
 if (y >= 0.d0) then
   a=acos(x/r)/(2.d0*pi)
 else
   a=1.d0-acos(x/r)/(2.d0*pi)
 end if

 end subroutine polarposition
!----------------------------!