!
!  This code integrates a two-dimensional linear dynamical system 
!  via predictor-corrector method (modified Euler)
!  This is the famous lovers' problem: Strogatz, H.S., Love affairs 
!  and differential equations, Math. Magazine, 61, 35 (1988).
!  The system has the general form
!  
!    dX      
!    -- = A X   
!    dt    
!
!  where 
!                    / a  b \               / x \
!              A =  |        |          X = |   | 
!                    \ c  d /               \ y /
!   
!
!  Instructions
!
!  - compile the code to produce an executable (.gfx extension)
!    
!          gfortran code.f90 -o code.gfx 
!  
!  - run the executable 
!
!          ./code.gfx 
!
!  - This will produce a three-column file (t,x,y) naamed
!    'lovers_destiny.dat'. Use gnuplot to plot your results. 
!    Type gnuplot at the prompt in a terminal and then type 
!    at the gnuplot prompt
!
!          p 'lovers_destiny.dat' u 2:3 w l 
!
!    to see the trajectory in the phase plane.
!    In gnuplot one uses abbreviations. For example
!
!    p   = plot
!    w l = with lines
!    pt 6 = pointtype 6 (open circles)
!    pt 7 = pointtype 7 (closed circles)
!
!    Refer to the website www.gnuplot.info for simple tutorials.
!
!================================================================
! Declarations
!
  integer, parameter :: kdp = kind(0.d0)
  real (kind=kdp) a,b,c,d
  complex  (kind=kdp) lam1,lam2,tr,det
  real (kind=kdp) dt, T_tot, x0, y0, xt, yt, x, y, ts, tk, L0, r
  integer N
!
! Parameters
!
!========================================================================
! A passionate and a fickle lover 

  b     =  0.8              !  Response of lover 1 to lover 2's feelings                 
  a     = 0.0               !  Response to one own's feelings (lover 1)
  c     = -b                !  Response of lover 2 to lover 1's feelings
  d     = 0.0               !  Response to one own's feelings (lover 2)
!========================================================================
! The two cautious lovers 

!  b     =  0.1             !  Response to the lover's feelings                 
!  a     = -0.3             !  Response to one own's feelings
!  c     = b
!  d     = a
!========================================================================  
  dt    = 0.001             ! Integration step
  dt2   = dt / 2.d0         ! half step
  T_tot = 10.               ! Total integration time (units of tau)
  N     = int ( T_tot/dt )  ! Number of integration steps
  ts    = 0.1               ! Sampling time (units of tau)
  
  tk    = 0.0               ! Initial value of running sampling time

  x0    = 0.1               ! Initial values
  y0    = 0.2               ! She loves him a little more at t = 0
                            
!========================================================================
! Compute spectrum of A
                   
  tr = a + d
  det = a*d - b*c
  lam1 = 0.5 * (tr + sqrt(tr**2 - 4*det) )
  lam2 = 0.5 * (tr - sqrt(tr**2 - 4*det) )    
  print *, ' lamda_1 =', lam1
  print *, ' lamda_2 =', lam2
!========================================================================
   
  open (1,file='lovers_destiny.dat')
  x = x0
  y = y0
  ! 
  !==================================================
  ! Modified Euler scheme : a simple Euler step
  ! is used to estimate the derivative as a centered 
  ! (midpoint) difference instead of using the 
  ! leftmost value
  !
  do i=1, N
     t = i*dt
     xt = x + (a*x + b*y) * dt
     yt = y + (c*x + d*y) * dt     
     x  = x + ( a*(x + xt) + b*(y + yt) )*dt2
     y  = y + ( c*(x + xt) + d*(y + yt) )*dt2     
     if ( t > tk ) then
         write (1,'(3(f12.8,3x))') t, x, y
         tk = tk + ts
     end if 
  end do
  !==================================================
   
  end program






