I need to perform an exercise with Fortran90 that has as its objective what is written in the title. I have a starting system formed by second degree equations with certain boundary conditions. I transformed the system to have equivalent first degree equations in order to integrate it, and then it was dimensionalized obtaining : dx/dt= vx; d(vx)/dt=-1; dy/dt= vy; d(vy)/dt=-1 with the following boundary conditions x(0)=1; y(0)=0; y(t)=0; x(t)=-1 with t=1.
Now I solve with the method of shooting and RK4 with the initial conditions so chosen: x(0)=1; y(0)=0; vx(0)=0; vy(0)=v1+Δv. Being a nonlinear system I have to iterate until convergence, the rendezvous is considered successfully completed at time t when the two vehicles are within 3 m, i.e. the set tolerance, to be dimensionalized. I report my code below:
SUBROUTINE dydx(neq, x, y, f)
INTEGER, INTENT(IN) :: neq
REAL*8, INTENT(IN) :: x
REAL*8, DIMENSION(neq), INTENT(IN) :: y !y=(y1,y2)=(1,v)
REAL*8, DIMENSION(neq),INTENT(OUT) :: f
!y 1 componente y, 2 componente v
f(1) = y(2) !dx/dt=v
f(2) = -y(1) !dv/dt=-1
END SUBROUTINE dydx
SUBROUTINE RK4(neq, h, x, yold, ynew)
INTEGER, INTENT(IN) :: neq
REAL*8, INTENT(IN) :: h, x
REAL*8, DIMENSION(neq), INTENT(IN) :: yold
REAL*8, DIMENSION(neq), INTENT(OUT) :: ynew
REAL*8, DIMENSION(neq) :: k1, k2, k3, k4
INTEGER :: i
CALL dydx(neq, x, yold, k1)
DO i=1, neq
ynew(i) = yold(i) + 0.5d0*h*k1(i)
END DO
CALL dydx(neq, x + 0.5d0*h, ynew, k2)
DO i=1, neq
ynew(i) = yold(i) + 0.5d0*h*k2(i)
END DO
CALL dydx(neq, x + 0.5d0*h, ynew, k3)
DO i=1, neq
ynew(i) = yold(i) + h*k3(i)
END DO
CALL dydx(neq, x + h, ynew, k4)
DO i=1, neq
ynew(i) = yold(i) + h*(k1(i) + 2.0d0*k2(i) + 2.0d0*k3(i) + k4(i)) / 6.0d0
END DO
END SUBROUTINE RK4
SUBROUTINE save_results(fname, neq, npoints, x, y)!!!!!risultati da salvare
CHARACTER (len=*), INTENT(IN) :: fname
INTEGER, INTENT(IN) :: npoints, neq
REAL*8, DIMENSION(0:npoints), INTENT(IN) :: x
REAL*8, DIMENSION(neq, 0:npoints), INTENT(IN) :: y
INTEGER :: i,j
OPEN(40, FILE=fname)
DO i=0, npoints
WRITE(40,'(5(1pe20.12))') x(i), (y(j,i), j=1, neq) !!!!!!
END DO
CLOSE(40)
END SUBROUTINE save_results
PROGRAM secondo
IMPLICIT NONE
INTEGER, PARAMETER :: npoints = 1000
INTEGER, PARAMETER :: neq = 4
INTEGER :: i, toll
REAL*8, DIMENSION(0:npoints) :: x
REAL*8, DIMENSION(neq,0:npoints) :: y
REAL*8 :: h, xmin, xmax
REAL*8 :: w1, w2, w3, y1, y2, yn, xn
!shooting
!condizioni al contorno
xmin = 0.0d0 !t0=0
xmax = 1.0d0 !t_rv=1
xn=-1.0d0
yn= 0.0d0
h = (xmax - xmin) / npoints
DO i=0, npoints
x(i) = xmin + i*h
END DO
!condizioni iniziali
y(1,0) = 1.0d0 !x(0)
y(2,0) = 0.0d0 !y(0)
y(3,0) = 0.0d0 !vx
y(4,0) = 0.5d0!vy
w1 = y(4,0)
DO i=1, npoints
CALL RK4(neq, h, x(i), y(:,i-1), y(:,i))
END DO
y1 = y(1,npoints)
CALL save_results('risultati_prima_iter.txt', neq, npoints, x, y)
y(1,0) = 1.0d0 !x(0)
y(2,0) = 0.0d0 !y(0)
y(3,0) = 0.0d0 !vx
y(4,0) = 2.0d0 !vy
w2 = y(4,0)
DO i=1, npoints
CALL RK4(neq, h, x(i), y(:,i-1), y(:,i))
END DO
y2 = y(1,npoints)
CALL save_results('risultati_seconda_iter.txt', neq, npoints, x, y)
w3 = w2 + (w2 -w1) / (y2 -y1) * (yn - y2)
!!!!!convergenza
toll=1.d-12
DO i=1, npoints
DO WHILE (ABS(yn-y2)<=toll)
CALL RK4(neq, h, x(i), y(:,i-1), y(:,i))
w2=w3
END DO
CALL save_results('risultati_terza_iter.txt', neq, npoints, x, y)
END DO
END PROGRAM secondo
At this point, if the code is correct, I should have found the value of Δv searched (in m/s) but I can not understand which column of the file "results_third_iter.txt" gives me the calculated result. In fact the next step is to make two graphs with gnuplot: the first with the trajectory of vehicle 1 in time, the second with the distance from the earth's surface of vehicle 1, as a function of time (distances in km, time in minutes) but even here I can not get sensible results.
Which data column in the file should I use to plot the graph? Also I think there is some error in the code because if I plot the three files at the same time without specifying the data columns to use, I get three identical parabolas which should not be so, where am I wrong?
Thanks to anyone who can help me
