# python script for generating data for Figure 1 in 
# Parametric resonance, chaos and spatial structure in the Lotka--Volterra model
# Mohamed Swailem and Alastair M. Rucklidge
# Copyright 2025 University of Leeds
# Unless otherwise stated, this dataset is licensed under a Creative Commons 
# Attribution 4.0 International Licence: 
# https://creativecommons.org/licenses/by/4.0/.

# solve Lotka-Volterra PDE with time-dependent carrying capacity
# uses ETD2 timestepping in spectral space (see Cox & Matthews J Comp Phys (2002) 176, 430)

import numpy as np 
import math 
import matplotlib.pyplot as plt 
import cmath 
from scipy.fftpack import fft, ifft 
from scipy.integrate import solve_ivp

# Version 1.00     # 17 April 2025: first working version to do multiple parameter values
#                  #   Didn't print version number
# version = 1.01   # added printing a version number
#version = 1.02    # added printing a summary line, increased initial mean field transient time to 100 periods
#                  # changed random amplitude to be relative rather than absolute
# version = 1.03   # added option of reading initial condition from a file, and running two periods only 5 May 2025
# version = 1.04   # changed stroboscopic map to 1 period, other tweaks
                   # changed stopping criterion to arms < 1.0e-12 * abs(abar)
                   # reduced the number of graphs produced
version = 1.05   # June 2025 -- added more control over initial condition for testing Turing instability
                 # added Ly, Da and Db to the file name
                 # added 1D as well as 2D
# edited for publication

pdf_png_extension = '.pdf'   # save as pdf or png
index = 0 # not used for a single calculation

# parameters: omega0, astar, kappa0, kappa1, n, alpha, Da, Db, domain size, time over which to solve

# adot = Da*laplacian(a) + lamda(t)*a*(b-1)                      # predator
# bdot = Db*laplacian(b) + b*(1-kappa(t)*(a+b)) - lamda(t)*a*b   # prey

# omega = omega0 / n
# lamda0    = (1 - (1+astar)*kappa0) / astar
# lamda1(t) = (1 - (1+astar)*(kappa0 + kappa1*cos(omega*t))) / astar
# lamda(t) = (1 - alpha)*lamda0 + alpha*lamda1 

# astar, kappa0, kappa1, n, alpha
kappa0 = 0.25
kappa1 = 0.24    # instability with alpha=1 is around kappa1=0.225
astar = 1
n = 0.6
alpha = 0.0

# derived quantities
lamda0 = (1 - (1 + astar) * kappa0) / astar
# lamda1 = (1 - (1 + astar) * (kappa0 + kappa1*cos(omega*time)) / astar
omega0squared = (1-kappa0)*(1-kappa0)/astar - kappa0*(1-kappa0) - kappa0*kappa0/4
omega0 = math.sqrt(omega0squared)
omega = omega0 / n

read_ic_from_file = 1
run_two_periods_only = 1

# compute a "mean-field" version at the same parameter values
def LV(t, ab):
    a, b = ab
    kappa = kappa0 + kappa1 * math.cos(omega*t)
    lamda1 = (1 - (1 + astar) * kappa) / astar
    lamda = (1-alpha) * lamda0 + alpha * lamda1
    if (lamda<0) or (kappa<0) or (omega0squared<0):
        print("Problem with parameters: lamda, kappa, omega0squared are:", lamda, kappa, omega0squared)
        error 
    return [lamda * a*(b - 1), b*(1 - (a + b) * kappa) - lamda * a * b]

# initial conditions
a_ic = astar * 1.05
b_ic = 1.05
random_amplitude = 1.0e-6

# run the mean field for one thousand periods to get over a transient, into the chaotic attractor
n_mean_field_periods_for_ic = 1000
print("Running the mean field ODE to generate initial condition: number of periods = ", n_mean_field_periods_for_ic)
sol = solve_ivp(LV, [0,n_mean_field_periods_for_ic * 2 * math.pi / omega], [a_ic, b_ic], method='RK45', rtol=1e-12, atol=1e-12, max_step=0.1)
mf_stroboscopic_a = sol.y[0]
mf_stroboscopic_b = sol.y[1]
a_ic = mf_stroboscopic_a[-1]
b_ic = mf_stroboscopic_b[-1]

# Da, Db  -- smaller than 0.01 causes trouble with the non-Taylor estimates of ETD2a etc
Da = 1.0
Db = 1.0

# domain size (we don't have a critical wavenumber a priori)
Lx = 500.0
Nx = 512
Ly = Lx
Ny = Nx  
# set Ly and Ny to zero for solving in 1D
# Lx = 1000.0
# Nx = 2048 * 4
# Ly = 0
# Ny = 0

# time stepping
timesteps_per_period = 100
timestep = 2*math.pi / omega / timesteps_per_period
timesteps_per_two_periods = 2 * timesteps_per_period
n_timesteps_plot = timesteps_per_period

if run_two_periods_only == 0:
    n_timesteps = 5000 * timesteps_per_two_periods
if run_two_periods_only == 1:
    n_timesteps = 1 * timesteps_per_two_periods
    n_timesteps_plot = timesteps_per_period

tmax = n_timesteps * timestep

# output control
save_figs_as_pdfs = 1

# set up arrays to store the summary output
save_time = np.zeros(n_timesteps+1)
save_abar = np.zeros(n_timesteps+1)
save_bbar = np.zeros(n_timesteps+1)
save_arms = np.zeros(n_timesteps+1)
save_brms = np.zeros(n_timesteps+1)
save_amax = np.zeros(n_timesteps+1)
save_amin = np.zeros(n_timesteps+1)
save_bmax = np.zeros(n_timesteps+1)
save_bmin = np.zeros(n_timesteps+1)
save_stroboscopic_time = np.zeros(n_timesteps//n_timesteps_plot + 1)
save_stroboscopic_abar = np.zeros(n_timesteps//n_timesteps_plot + 1)
save_stroboscopic_bbar = np.zeros(n_timesteps//n_timesteps_plot + 1)
save_stroboscopic_arms = np.zeros(n_timesteps//n_timesteps_plot + 1)
save_stroboscopic_brms = np.zeros(n_timesteps//n_timesteps_plot + 1)
save_stroboscopic_amax = np.zeros(n_timesteps//n_timesteps_plot + 1)
save_stroboscopic_amin = np.zeros(n_timesteps//n_timesteps_plot + 1)
save_stroboscopic_bmax = np.zeros(n_timesteps//n_timesteps_plot + 1)
save_stroboscopic_bmin = np.zeros(n_timesteps//n_timesteps_plot + 1)
save_stroboscopic_wave = np.zeros(n_timesteps//n_timesteps_plot + 1)

# say what we are going to do
print("astar, kappa0, kappa1, n, alpha:", astar, kappa0, kappa1, n, alpha)
print("lamda0, omega0, omega:", lamda0, omega0, omega)
print("a_ic, b_ic, random_amplitude:", a_ic, b_ic, random_amplitude)
print("timesteps_per_period, n_timesteps:", timesteps_per_period, n_timesteps)
print("Lx, Ly, Nx, Ny, Da, Db:", Lx, Ly, Nx, Ny, Da, Db)

# construct a filename
file_name_1 = "LV_PDE_"
file_name_2 = "alpha_" + "{:06.4f}".format(alpha) + \
              "_kappa0_" + "{:06.4f}".format(kappa0) + \
              "_kappa1_" + "{:06.4f}".format(kappa1) + \
              "_n_" + "{:06.4f}".format(n) + \
              "_Da_" + "{:08.4f}".format(Da) + \
              "_Db_" + "{:08.4f}".format(Db) + \
              "_Lx_" + "{:06.1f}".format(Lx) + \
              "_Ly_" + "{:06.1f}".format(Ly) + \
              "_" + "{:04d}".format(index)              
file_name_2 = file_name_2.replace(".", "p") 
print("File name: ", file_name_2)

# setup the (x,y) and (kx,ky) coordinates, 1D first
x=np.arange(0.0,Nx)*(Lx/Nx) 
if Ny != 0:
    y=np.arange(0.0,Ny)*(Ly/Ny) 

kx=np.arange(0,Nx) * 1.0
if Ny != 0:
    ky=np.arange(0,Ny) * 1.0

for i in range(int(Nx/2+1),Nx): 
    kx[i]= i-Nx 

if Ny != 0:
    for j in range(int(Ny/2+1),Ny): 
        ky[j]=j-Ny 

Kx=kx*((2*math.pi)/Lx) 
if Ny != 0:
    Ky=ky*((2*math.pi)/Ly) 

# and now a 2D grid of them (could use meshgrid for this)
if Ny != 0:
    xx=np.zeros((Nx,Ny)) 
    yy=np.zeros((Nx,Ny)) 
    KKx=np.zeros((Nx,Ny)) 
    KKy=np.zeros((Nx,Ny)) 
else:
    xx=np.zeros((Nx)) 
    KKx=np.zeros((Nx)) 

if Ny != 0:
    for i in range(0,Nx): 
        for j in range(0,Ny): 
            xx[i,j]=x[i] 
            yy[i,j]=y[j] 
            KKx[i,j]=Kx[i] 
            KKy[i,j]=Ky[j] 
else:
    for i in range(0,Nx): 
        xx[i]=x[i] 
        KKx[i]=Kx[i] 

if Ny != 0:
    KK2=pow(KKx,2)+pow(KKy,2)            # k^2 for the Laplacian
else:
    KK2=pow(KKx,2)                       # k^2 for the Laplacian

# now set up the arrays needed for the ETD2 method, 
# from Cox and Matthews, J Comp Phys 2002

L_a = - Da * KK2          # the linear part of the PDE
L_b = - Db * KK2
h = timestep
exp_hL_a = np.exp(h * L_a) 
exp_hL_b = np.exp(h * L_b) 

# for these, don't divide by c^2h until the end 
ETD2a_a = (exp_hL_a*(1 + L_a*h) - 1 - 2*L_a*h )   #  / L_a / L_a / h
ETD2a_b = (exp_hL_b*(1 + L_b*h) - 1 - 2*L_b*h )   #  / L_b / L_b / h
ETD2b_a = (exp_hL_a*(  -1) + 1 + L_a*h )          #  / L_a / L_a / h
ETD2b_b = (exp_hL_b*(  -1) + 1 + L_b*h )          #  / L_b / L_b / h

# the method here is:
# Uhat_next = exp_hL * Uhat + ETD2a * NLhat + ETD2b * NLoldhat
# NB there will be an error if L is zero (or small)
# in this case, alternate methods should be used to compute these
# Taylor series in the case where L=0 gives ETD=h, ETD2a=3*h/2, ETD2b=-h/2
# properly, we should include a small-L expansion here
if Ny != 0:
    for i in range(0,Nx): 
        for j in range(0,Ny): 
                if (L_a[i,j] == 0.0) :
                        ETD2a_a[i,j] = (3*h)/2
                        ETD2b_a[i,j] = -h/2
                else:
                        ETD2a_a[i,j] = ETD2a_a[i,j] / L_a[i,j] / L_a[i,j] / h
                        ETD2b_a[i,j] = ETD2b_a[i,j] / L_a[i,j] / L_a[i,j] / h
else:
    for i in range(0,Nx): 
        if (L_a[i] == 0.0) :
            ETD2a_a[i] = (3*h)/2
            ETD2b_a[i] = -h/2
        else:
            ETD2a_a[i] = ETD2a_a[i] / L_a[i] / L_a[i] / h
            ETD2b_a[i] = ETD2b_a[i] / L_a[i] / L_a[i] / h
 
if Ny != 0:
    for i in range(0,Nx): 
        for j in range(0,Ny): 
                if (L_b[i,j] == 0.0) :
                        ETD2a_b[i,j] = (3*h)/2
                        ETD2b_b[i,j] = -h/2
                else:
                        ETD2a_b[i,j] = ETD2a_b[i,j] / L_b[i,j] / L_b[i,j] / h
                        ETD2b_b[i,j] = ETD2b_b[i,j] / L_b[i,j] / L_b[i,j] / h
else:
    for i in range(0,Nx): 
        if (L_b[i] == 0.0) :
            ETD2a_b[i] = (3*h)/2
            ETD2b_b[i] = -h/2
        else:
            ETD2a_b[i] = ETD2a_b[i] / L_b[i] / L_b[i] / h
            ETD2b_b[i] = ETD2b_b[i] / L_b[i] / L_b[i] / h

# initial conditions
if read_ic_from_file == 0:
    print("Generating random initial conditions, after transient in mean field ODEs")
    if Ny != 0:
        ap = a_ic * (1 + random_amplitude*np.random.rand(Nx,Ny))
        bp = b_ic * (1 + random_amplitude*np.random.rand(Nx,Ny))
    else:
        ap = a_ic * (1 + random_amplitude*np.random.rand(Nx))
        bp = b_ic * (1 + random_amplitude*np.random.rand(Nx))
if read_ic_from_file == 1:
    print("Reading initial conditions from: ", file_name_1 + file_name_2 + "_a_ic.txt.gz")
    ap = np.loadtxt(file_name_1 + file_name_2 + "_a_ic.txt.gz")
    bp = np.loadtxt(file_name_1 + file_name_2 + "_b_ic.txt.gz")

# we will work in spectral space for the timestepping
ahat = np.fft.fftn(ap) 
bhat = np.fft.fftn(bp) 

# smooth the random initial conditions a little, so we are not differentiating noise
if read_ic_from_file == 0:
    if Ny != 0:
        for i in range(0,Nx): 
            for j in range(0,Ny): 
                ahat[i,j]=ahat[i,j] / (1+0.01*KK2[i,j]*KK2[i,j])
                bhat[i,j]=bhat[i,j] / (1+0.01*KK2[i,j]*KK2[i,j])
    else:
        for i in range(0,Nx): 
            ahat[i]=ahat[i] / (1+0.01*KK2[i]*KK2[i])
            bhat[i]=bhat[i] / (1+0.01*KK2[i]*KK2[i])

ap = np.real(np.fft.ifftn(ahat))
bp = np.real(np.fft.ifftn(bhat))

# # display the initial condition in physical space -- only if run_two_periods_only
if run_two_periods_only == 1:
    fig_contour, fig_contour_axes = plt.subplots(nrows = 1, ncols = 1, figsize=(8, 6), tight_layout=True)
    if Ny != 0:
        fig_contour_contour = fig_contour_axes.contour(xx, yy, ap) 
        fig_contour_axes.set_aspect('equal')
        fig_contour_axes.set_title("Contour of initial a(x,y)")
        fig_contour_axes.set_xlabel('$x$')
        fig_contour_axes.set_ylabel('$y$')
        fig_contour.colorbar(fig_contour_contour, location = 'right')
    else:
        fig_contour_axes.plot(x, ap)
        fig_contour_axes.set_title("Plot of initial a(x)")
        fig_contour_axes.set_xlabel('$x$')
        fig_contour_axes.set_ylabel('$ap$')
    fig_contour.tight_layout()

    if save_figs_as_pdfs == 1:
        fig_contour.savefig(file_name_1 + file_name_2 + "_contour_a_IC" + pdf_png_extension, bbox_inches='tight') 
    else:
        plt.show(Block=True)
    plt.close(fig_contour)

# and now the timestep you've all be waiting for...
# we are timestepping in spectral space,
# so we start with ahat and bhat from the previous timestep

# adot = Da*laplacian(a) + lamda(t)*a*(b-1)                      # predator
# bdot = Db*laplacian(b) + b*(1-kappa(t)*(a+b)) - lamda(t)*a*b   # prey
# lamda0    = (1 - (1+astar)*kappa0) / astar
# lamda1(t) = (1 - (1+astar)*(kappa0 + kappa1*cos(omega*t))) / astar
# lamda(t) = (1 - alpha)*lamda0 + alpha*lamda1 

def Diagnostics(ap, bp):
    if Ny != 0:
        abar = np.sum(ap) / Nx / Ny
        bbar = np.sum(bp) / Nx / Ny
        arms = math.sqrt(np.sum(np.square(ap - abar)) / Nx / Ny)
        brms = math.sqrt(np.sum(np.square(bp - bbar)) / Nx / Ny)
    else:
        abar = np.sum(ap) / Nx
        bbar = np.sum(bp) / Nx
        arms = math.sqrt(np.sum(np.square(ap - abar)) / Nx)
        brms = math.sqrt(np.sum(np.square(bp - bbar)) / Nx)
    amax = np.max(ap)
    amin = np.min(ap)
    bmax = np.max(bp)
    bmin = np.min(bp)
    return [abar, bbar, arms, brms, amax, amin, bmax, bmin]

stroboscopic_step = 0
for istep in range(0, n_timesteps): 
        the_time = istep * timestep  # this is the time at the start of the timestep
        kappa = kappa0 + kappa1*math.cos(omega*the_time)
        lamda1 = (1 - (1+astar)*kappa) / astar
        lamda = (1 - alpha)*lamda0 + alpha*lamda1

        ap = np.real(np.fft.ifftn(ahat))
        bp = np.real(np.fft.ifftn(bhat))

# Diagnostics -- defined above -- we have the fields in physical space at the start of each timestep.
        [abar, bbar, arms, brms, amax, amin, bmax, bmin] = Diagnostics(ap, bp)

# save and print out our progress
        save_time[istep] = the_time
        save_abar[istep] = abar
        save_bbar[istep] = bbar
        save_arms[istep] = arms
        save_brms[istep] = brms
        save_amax[istep] = amax
        save_amin[istep] = amin
        save_bmax[istep] = bmax
        save_bmin[istep] = bmin

        stop_after_this_step = 0

        if (istep % n_timesteps_plot == 0): 
                penultimate_ap = ap # save this each time so we have the previous one to plot as well
                save_stroboscopic_time[stroboscopic_step] = the_time
                save_stroboscopic_abar[stroboscopic_step] = abar    
                save_stroboscopic_bbar[stroboscopic_step] = bbar    
                save_stroboscopic_arms[stroboscopic_step] = arms    
                save_stroboscopic_brms[stroboscopic_step] = brms
                save_stroboscopic_amax[stroboscopic_step] = amax
                save_stroboscopic_amin[stroboscopic_step] = amin
                save_stroboscopic_bmax[stroboscopic_step] = bmax
                save_stroboscopic_bmin[stroboscopic_step] = bmin
# estimate a lengthscale from the Laplacian: From Bentley et al: local wavenumber = sqrt(-Lapap/a)
                Lapap = np.real(np.fft.ifftn( - KK2 * ahat)) 
                local_wavenumber_squared = np.abs(- Lapap / ap)
                if Ny != 0:
                    local_wavenumber_rms = math.sqrt(np.sum(local_wavenumber_squared) / Nx / Ny)
                else:
                    local_wavenumber_rms = math.sqrt(np.sum(local_wavenumber_squared) / Nx)
                save_stroboscopic_wave[stroboscopic_step] = local_wavenumber_rms
                print("{:4d}".format(stroboscopic_step), "{:6d}".format(istep), "{:8.1f}".format(the_time * omega/2/math.pi), "{:20.6f}".format(abar), "{:20.6g}".format(arms), "{:20.6g}".format(local_wavenumber_rms))
                stroboscopic_step = stroboscopic_step + 1
                if arms < 1.0e-12 * abs(abar):          # commenting this out can lead to errors
                    stop_after_this_step = istep        # in the contour plot

# Now do the timestep: compute nonlinear terms in physical space and in Fourier space
# ODE:    return [lamda * a*(b - 1), b*(1 - (a + b) * kappa) - lamda * a * b]
        NL_a = lamda * ap * (bp - 1)
        NL_b = bp * (1 - kappa * (ap + bp)) - lamda * ap * bp

        NL_ahat = np.fft.fftn(NL_a)
        NL_bhat = np.fft.fftn(NL_b)
# first step is special: setting NLoldhat=NLhat is equivalent to ETD1
        if (istep == 0):
                NL_aoldhat = NL_ahat
                NL_boldhat = NL_bhat

        ahat_next = exp_hL_a * ahat + ETD2a_a * NL_ahat + ETD2b_a * NL_aoldhat
        bhat_next = exp_hL_b * bhat + ETD2a_b * NL_bhat + ETD2b_b * NL_boldhat
        if stop_after_this_step == 0:
                                  # only update these if we are progressing
                                  # if we want to stop, we want ahat, bhat at the start of this timestep
            ahat = ahat_next      # whereas these are now at the end of the timestep
            bhat = bhat_next
            NL_aoldhat = NL_ahat
            NL_boldhat = NL_bhat

        istep_save = istep
        if stop_after_this_step != 0:
                break

if stop_after_this_step != 0:
    print("Stopping at timestep ", istep_save)
else:
    istep_save = istep_save + 1

istep = istep_save
the_time = istep * timestep

ap = np.real(np.fft.ifftn(ahat))
bp = np.real(np.fft.ifftn(bhat))

[abar, bbar, arms, brms, amax, amin, bmax, bmin] = Diagnostics(ap, bp)

# save and print out our progress
save_time[istep] = the_time
save_abar[istep] = abar
save_bbar[istep] = bbar
save_arms[istep] = arms
save_brms[istep] = brms
save_amax[istep] = amax
save_amin[istep] = amin
save_bmax[istep] = bmax
save_bmin[istep] = bmin

save_stroboscopic_time[stroboscopic_step] = the_time
save_stroboscopic_abar[stroboscopic_step] = abar    
save_stroboscopic_bbar[stroboscopic_step] = bbar    
save_stroboscopic_arms[stroboscopic_step] = arms    
save_stroboscopic_brms[stroboscopic_step] = brms    
save_stroboscopic_amax[stroboscopic_step] = amax
save_stroboscopic_amin[stroboscopic_step] = amin
save_stroboscopic_bmax[stroboscopic_step] = bmax
save_stroboscopic_bmin[stroboscopic_step] = bmin

# estimate a lengthscale from the Laplacian: From Bentley et al: local wavenumber = sqrt(-Lapap/a)
Lapap = np.real(np.fft.ifftn( - KK2 * ahat)) 
local_wavenumber_squared = np.abs(- Lapap / ap)
if Ny != 0:
    local_wavenumber_rms = math.sqrt(np.sum(local_wavenumber_squared) / Nx / Ny)
else:
    local_wavenumber_rms = math.sqrt(np.sum(local_wavenumber_squared) / Nx)
save_stroboscopic_wave[stroboscopic_step] = local_wavenumber_rms
print("{:4d}".format(stroboscopic_step), "{:6d}".format(istep), "{:8.1f}".format(the_time * omega/2/math.pi), "{:20.6f}".format(abar), "{:20.6g}".format(arms), "{:20.6g}".format(local_wavenumber_rms))
stroboscopic_step = stroboscopic_step + 1

# compute a "mean-field" version at the same parameter values
t_eval = 2 * math.pi / omega * np.arange(stroboscopic_step) * 2 
sol = solve_ivp(LV, [0,t_eval[-1]], [a_ic, b_ic], t_eval=t_eval, method='RK45', rtol=1e-12, atol=1e-12, max_step=0.1)
        
if sol.status != 0:
    print("Status: ", sol.status)
        
mf_stroboscopic_time = sol.t
mf_stroboscopic_a = sol.y[0]
mf_stroboscopic_b = sol.y[1]

# save / display the final state in physical space (but only if we didn't converge to a flat state)
if stop_after_this_step == 0:

    np.savetxt(file_name_1 + file_name_2 + "_a.txt.gz", ap)
    np.savetxt(file_name_1 + file_name_2 + "_b.txt.gz", bp)
# save penultimate ap as well, but only if running for two periods
    if run_two_periods_only == 1:
        np.savetxt(file_name_1 + file_name_2 + "_penultimate_a.txt.gz", penultimate_ap)

    fig_contour, fig_contour_axes = plt.subplots(nrows = 1, ncols = 1, figsize=(8, 6), tight_layout=True)
    if Ny != 0:
        fig_contour_contour = fig_contour_axes.contour(xx, yy, ap) 
        fig_contour_axes.set_aspect('equal')
        fig_contour_axes.set_title("Contour of final a(x,y)")
        fig_contour_axes.set_xlabel('$x$')
        fig_contour_axes.set_ylabel('$y$')
        fig_contour.colorbar(fig_contour_contour, location = 'right')
    else:
        fig_contour_axes.plot(x, ap)
        fig_contour_axes.set_title("Plot of final a(x)")
        fig_contour_axes.set_xlabel('$x$')
        fig_contour_axes.set_ylabel('$ap$')

    fig_contour.tight_layout()
    if save_figs_as_pdfs == 1:
        fig_contour.savefig(file_name_1 + file_name_2 + "_contour_a_final" + pdf_png_extension, bbox_inches='tight') 
    else:
        plt.show(Block=True)
    plt.close(fig_contour)

# display the stroboscopic phase portrait in log coordinates
    fig_phase_portrait_stroboscopic, fig_phase_portrait_stroboscopic_axes = plt.subplots(nrows = 1, ncols = 1, figsize=(8, 6), tight_layout=True)
    fig_phase_portrait_stroboscopic_axes.plot(np.log10(save_stroboscopic_abar[0:stroboscopic_step-1]), np.log10(save_stroboscopic_bbar[0:stroboscopic_step-1]), 'm+') 
    fig_phase_portrait_stroboscopic_axes.plot(np.log10(mf_stroboscopic_a), np.log10(mf_stroboscopic_b), 'c+') 
    fig_phase_portrait_stroboscopic_axes.set_title("phase portrait (stroboscopic) of (abar, bbar) plus mean field")
    fig_phase_portrait_stroboscopic_axes.set_xlabel('$log10 abar$')
    fig_phase_portrait_stroboscopic_axes.set_ylabel('$log10 bbar$')
    fig_phase_portrait_stroboscopic.tight_layout()
    if save_figs_as_pdfs == 1:
            fig_phase_portrait_stroboscopic.savefig(file_name_1 + file_name_2 + "_phase_portrait_stroboscopic" + pdf_png_extension, bbox_inches='tight') 
    else:
            plt.show(Block=True)
    plt.close(fig_phase_portrait_stroboscopic)

# display penultimate ap as well, but only if running for two periods
    if run_two_periods_only == 1:
        fig_contour, fig_contour_axes = plt.subplots(nrows = 1, ncols = 1, figsize=(8, 6), tight_layout=True)
        if Ny != 0:
            fig_contour_contour = fig_contour_axes.contour(xx, yy, penultimate_ap) 
            fig_contour_axes.set_aspect('equal')
            fig_contour_axes.set_title("Contour of penultimate a(x,y)")
            fig_contour_axes.set_xlabel('$x$')
            fig_contour_axes.set_ylabel('$y$')
            fig_contour.colorbar(fig_contour_contour, location = 'right')
        else:
            fig_contour_axes.plot(x, ap)
            fig_contour_axes.set_title("Plot of penultimate a(x)")
            fig_contour_axes.set_xlabel('$x$')
            fig_contour_axes.set_ylabel('$ap$')
        fig_contour.tight_layout()
        if save_figs_as_pdfs == 1:
            fig_contour.savefig(file_name_1 + file_name_2 + "_contour_a_penultimate" + pdf_png_extension, bbox_inches='tight') 
        else:
            plt.show(Block=True)
        plt.close(fig_contour)

#    fig_contour, fig_contour_axes = plt.subplots(nrows = 1, ncols = 1, figsize=(8, 6), tight_layout=True)
#    fig_contour_contour = fig_contour_axes.contour(xx, yy, bp) 
#    fig_contour_axes.set_aspect('equal')
#    fig_contour_axes.set_title("Contour of final b(x,y)")
#    fig_contour_axes.set_xlabel('$x$')
#    fig_contour_axes.set_ylabel('$y$')
#    fig_contour.tight_layout()
#    fig_contour.colorbar(fig_contour_contour, location = 'right')
#    if save_figs_as_pdfs == 1:
#        fig_contour.savefig(file_name_1 + file_name_2 + "_contour_b_final" + pdf_png_extension, bbox_inches='tight') 
#    else:
#        plt.show(Block=True)
#    plt.close(fig_contour)

# display the final state in Fourier space
#    fig_fourier_spectrum, fig_fourier_spectrum_axes = plt.subplots(nrows = 1, ncols = 1, figsize=(8, 6), tight_layout=True)
#    ahat0 = ahat
#    ahat[0,0] = 1.0
#    fig_fourier_spectrum_axes.plot(np.sqrt(KK2), np.log10(np.abs(ahat0)), 'k+') 
#    fig_fourier_spectrum_axes.set_title("Fourier spectrum of final a(x,y)")
#    fig_fourier_spectrum_axes.set_xlabel('$k$')
#    fig_fourier_spectrum_axes.set_ylabel('$log10|ahat|$')
#    fig_fourier_spectrum.tight_layout()
#    if save_figs_as_pdfs == 1:
#        fig_fourier_spectrum.savefig(file_name_1 + file_name_2 + "_Fourier_spectrum_a_final" + pdf_png_extension, bbox_inches='tight') 
#    else:
#        plt.show(Block=True)
#    plt.close(fig_fourier_spectrum)

the_time             = save_stroboscopic_time[stroboscopic_step-1] 
abar                 = save_stroboscopic_abar[stroboscopic_step-1] 
bbar                 = save_stroboscopic_bbar[stroboscopic_step-1] 
arms                 = save_stroboscopic_arms[stroboscopic_step-1] 
brms                 = save_stroboscopic_brms[stroboscopic_step-1] 
amax                 = save_stroboscopic_amax[stroboscopic_step-1] 
amin                 = save_stroboscopic_amin[stroboscopic_step-1] 
bmax                 = save_stroboscopic_bmax[stroboscopic_step-1] 
bmin                 = save_stroboscopic_bmin[stroboscopic_step-1] 
local_wavenumber_rms = save_stroboscopic_wave[stroboscopic_step-1] 

print("Summary: ", file_name_2, " ", \
      "{:06.4f}".format(astar), "{:06.4f}".format(kappa0), "{:06.4f}".format(kappa1), "{:06.4f}".format(n), "{:06.4f}".format(alpha), \
      "{:06.4f}".format(lamda0), "{:06.4f}".format(omega0), "{:06.4f}".format(omega), \
      "{:06.4f}".format(a_ic), "{:06.4f}".format(b_ic), "{:06.4f}".format(random_amplitude), \
      "{:6d}".format(timesteps_per_period), "{:10d}".format(n_timesteps), \
      "{:06.4f}".format(Lx), "{:06.4f}".format(Ly), "{:6d}".format(Nx), "{:6d}".format(Ny), "{:06.4f}".format(Da), "{:06.4f}".format(Db), \
      "{:5.2f}".format(version), \
      "{:20.6g}".format(the_time * omega / 2/math.pi), \
      "{:20.6g}".format(abar), "{:20.6g}".format(bbar), "{:20.6g}".format(arms), "{:20.6g}".format(brms), \
      "{:20.6g}".format(amax), "{:20.6g}".format(amin), "{:20.6g}".format(bmax), "{:20.6g}".format(bmin), \
      "{:20.6g}".format(local_wavenumber_rms))

# display the phase portrait
fig_phase_portrait, fig_phase_portrait_axes = plt.subplots(nrows = 1, ncols = 1, figsize=(8, 6), tight_layout=True)
fig_phase_portrait_axes.plot(np.log10(save_abar[0:istep_save]), np.log10(save_bbar[0:istep_save]), '0.8') 
fig_phase_portrait_axes.plot(np.log10(save_stroboscopic_abar[0:stroboscopic_step-1]), np.log10(save_stroboscopic_bbar[0:stroboscopic_step-1]), 'm+') 
fig_phase_portrait_axes.set_title("phase_portrait of (abar, bbar)")
fig_phase_portrait_axes.set_xlabel('$log10 abar$')
fig_phase_portrait_axes.set_ylabel('$log10 bbar$')
fig_phase_portrait.tight_layout()
if save_figs_as_pdfs == 1:
        fig_phase_portrait.savefig(file_name_1 + file_name_2 + "_phase_portrait" + pdf_png_extension, bbox_inches='tight') 
else:
        plt.show(Block=True)
plt.close(fig_phase_portrait)

# display the time dependence of the RMS variation in log coordinates
fig_time_series_RMS_variation, fig_time_series_RMS_variation_axes = plt.subplots(nrows = 1, ncols = 1, figsize=(12, 6), tight_layout=True)
fig_time_series_RMS_variation_axes.plot(save_time[0:istep_save] * omega/2/math.pi, np.log10(save_arms[0:istep_save]), '0.8') 
fig_time_series_RMS_variation_axes.plot(save_stroboscopic_time[0:stroboscopic_step-1] * omega/2/math.pi, np.log10(save_stroboscopic_arms[0:stroboscopic_step-1]), 'm+') 
fig_time_series_RMS_variation_axes.plot(save_stroboscopic_time[0:stroboscopic_step-1] * omega/2/math.pi, np.log10(save_stroboscopic_wave[0:stroboscopic_step-1]), 'g+') 
fig_time_series_RMS_variation_axes.set_title("Time dependence of the RMS variation of a")
fig_time_series_RMS_variation_axes.set_xlabel('$time/period$')
fig_time_series_RMS_variation_axes.set_ylabel('$log10 arms$')
fig_time_series_RMS_variation.tight_layout()
if save_figs_as_pdfs == 1:
        fig_time_series_RMS_variation.savefig(file_name_1 + file_name_2 + "_time_series_RMS_variation" + pdf_png_extension, bbox_inches='tight') 
else:
        plt.show(Block=True)
plt.close(fig_time_series_RMS_variation)

# # display the time dependence of the max/avg/min variation in log coordinates
# fig_time_series_max_bar_min_variation, fig_time_series_max_bar_min_variation_axes = plt.subplots(nrows = 1, ncols = 1, figsize=(12, 6), tight_layout=True)
# fig_time_series_max_bar_min_variation_axes.plot(save_time[0:istep_save] * omega/2/math.pi, np.log10(save_amax[0:istep_save]), '0.8') 
# fig_time_series_max_bar_min_variation_axes.plot(save_time[0:istep_save] * omega/2/math.pi, np.log10(save_abar[0:istep_save]), 'k') 
# fig_time_series_max_bar_min_variation_axes.plot(save_time[0:istep_save] * omega/2/math.pi, np.log10(save_amin[0:istep_save]), '0.8') 
# fig_time_series_max_bar_min_variation_axes.plot(save_stroboscopic_time[0:stroboscopic_step-1] * omega/2/math.pi, np.log10(save_stroboscopic_abar[0:stroboscopic_step-1]), 'm+') 
# fig_time_series_max_bar_min_variation_axes.set_title("Time dependence of the max/avg/min variation of a")
# fig_time_series_max_bar_min_variation_axes.set_xlabel('$time/period$')
# fig_time_series_max_bar_min_variation_axes.set_ylabel('$log10 max/avg/min$')
# fig_time_series_max_bar_min_variation.tight_layout()
# if save_figs_as_pdfs == 1:
#         fig_time_series_max_bar_min_variation.savefig(file_name_1 + file_name_2 + "_time_series_max_bar_min_variation" + pdf_png_extension, bbox_inches='tight') 
# else:
#         plt.show(Block=True)
# plt.close(fig_time_series_max_bar_min_variation)

# save all the data in files
LV_PDE_data_file = open(file_name_1 + file_name_2 + "_data.txt", mode = "wt")
print("# astar, kappa0, kappa1, n, alpha:", astar, kappa0, kappa1, n, alpha, file = LV_PDE_data_file)
print("# lamda0, omega0, omega:", lamda0, omega0, omega, file = LV_PDE_data_file)
print("# a_ic, b_ic, random_amplitude:", a_ic, b_ic, random_amplitude, file = LV_PDE_data_file)
print("# timesteps_per_period, n_timesteps:", timesteps_per_period, n_timesteps, file = LV_PDE_data_file)
print("# Lx, Ly, Nx, Ny, Da, Db:", Lx, Ly, Nx, Ny, Da, Db, file = LV_PDE_data_file)
print("# Version:", "{:5.2f}".format(version), file = LV_PDE_data_file)
print("#", file = LV_PDE_data_file)
print("#strob    istep          time/period                 abar                 bbar                 arms                 brms                 amax                 bmax                 amin                 bmin           wavenumber", file = LV_PDE_data_file)

stroboscopic_step = 0
for istep in range(0, istep_save + 1): 
        if (istep % n_timesteps_plot != 0): 
                print("{:6d}".format(-1), \
                      "{:8d}".format(istep), \
                      "{:20.6f}".format(save_time[istep] * omega/2/math.pi), \
                      "{:20.6f}".format(save_abar[istep]), \
                      "{:20.6f}".format(save_bbar[istep]), \
                      "{:20.6g}".format(save_arms[istep]), \
                      "{:20.6g}".format(save_brms[istep]), \
                      "{:20.6g}".format(save_amax[istep]), \
                      "{:20.6g}".format(save_bmax[istep]), \
                      "{:20.6g}".format(save_amin[istep]), \
                      "{:20.6g}".format(save_bmin[istep]), \
                      "{:20.6g}".format(-1), \
                      file = LV_PDE_data_file)
        if (istep % n_timesteps_plot == 0): 
                print("{:6d}".format(stroboscopic_step), \
                      "{:8d}".format(istep), \
                      "{:20.6f}".format(save_time[istep] * omega/2/math.pi), \
                      "{:20.6f}".format(save_abar[istep]), \
                      "{:20.6f}".format(save_bbar[istep]), \
                      "{:20.6g}".format(save_arms[istep]), \
                      "{:20.6g}".format(save_brms[istep]), \
                      "{:20.6g}".format(save_amax[istep]), \
                      "{:20.6g}".format(save_bmax[istep]), \
                      "{:20.6g}".format(save_amin[istep]), \
                      "{:20.6g}".format(save_bmin[istep]), \
                      "{:20.6g}".format(save_stroboscopic_wave[stroboscopic_step]), \
                      file = LV_PDE_data_file)
                stroboscopic_step = stroboscopic_step + 1
LV_PDE_data_file.close()

# save all the data in files -- this one only save stroboscopic data to save space
LV_PDE_stroboscopic_data_file = open(file_name_1 + file_name_2 + "_stroboscopic_data.txt", mode = "wt")
print("# astar, kappa0, kappa1, n, alpha:", astar, kappa0, kappa1, n, alpha, file = LV_PDE_stroboscopic_data_file)
print("# lamda0, omega0, omega:", lamda0, omega0, omega, file = LV_PDE_stroboscopic_data_file)
print("# a_ic, b_ic, random_amplitude:", a_ic, b_ic, random_amplitude, file = LV_PDE_stroboscopic_data_file)
print("# timesteps_per_period, n_timesteps:", timesteps_per_period, n_timesteps, file = LV_PDE_stroboscopic_data_file)
print("# Lx, Ly, Nx, Ny, Da, Db:", Lx, Ly, Nx, Ny, Da, Db, file = LV_PDE_stroboscopic_data_file)
print("# Version:", "{:5.2f}".format(version), file = LV_PDE_stroboscopic_data_file)
print("#", file = LV_PDE_stroboscopic_data_file)
print("#strob    istep          time/period                 abar                 bbar                 arms                 brms                 amax                 bmax                 amin                 bmin           wavenumber", file = LV_PDE_stroboscopic_data_file)

stroboscopic_step = 0
for istep in range(0, istep_save + 1): 
        if (istep % n_timesteps_plot == 0): 
                print("{:6d}".format(stroboscopic_step), \
                      "{:8d}".format(istep), \
                      "{:20.6f}".format(save_time[istep] * omega/2/math.pi), \
                      "{:20.6f}".format(save_abar[istep]), \
                      "{:20.6f}".format(save_bbar[istep]), \
                      "{:20.6g}".format(save_arms[istep]), \
                      "{:20.6g}".format(save_brms[istep]), \
                      "{:20.6g}".format(save_amax[istep]), \
                      "{:20.6g}".format(save_bmax[istep]), \
                      "{:20.6g}".format(save_amin[istep]), \
                      "{:20.6g}".format(save_bmin[istep]), \
                      "{:20.6g}".format(save_stroboscopic_wave[stroboscopic_step]), \
                      file = LV_PDE_stroboscopic_data_file)
                stroboscopic_step = stroboscopic_step + 1
LV_PDE_stroboscopic_data_file.close()

