# 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/.

# reread the a data and replot the contour diagrams

import numpy as np 
import math 
import matplotlib.pyplot as plt 
from matplotlib.ticker import FormatStrFormatter

version = 1.00
# edited for publication

filename_n0p6_0 = 'LV_PDE_alpha_0p0000_kappa0_0p2500_kappa1_0p2400_n_0p6000_Da_001p0000_Db_001p0000_Lx_0500p0_Ly_0500p0_0000_penultimate_a'
filename_n0p6_1 = 'LV_PDE_alpha_0p0000_kappa0_0p2500_kappa1_0p2400_n_0p6000_Da_001p0000_Db_001p0000_Lx_0500p0_Ly_0500p0_0000_a'
filename_n0p7_0 = 'LV_PDE_alpha_0p0000_kappa0_0p2500_kappa1_0p2400_n_0p7000_Da_001p0000_Db_001p0000_Lx_0500p0_Ly_0500p0_0000_penultimate_a'
filename_n0p7_1 = 'LV_PDE_alpha_0p0000_kappa0_0p2500_kappa1_0p2400_n_0p7000_Da_001p0000_Db_001p0000_Lx_0500p0_Ly_0500p0_0000_a'

ap_n0p6_0 = np.loadtxt(filename_n0p6_0 + '.txt.gz')
ap_n0p6_1 = np.loadtxt(filename_n0p6_1 + '.txt.gz')
ap_n0p7_0 = np.loadtxt(filename_n0p7_0 + '.txt.gz')
ap_n0p7_1 = np.loadtxt(filename_n0p7_1 + '.txt.gz')

filename_list = [filename_n0p6_0, filename_n0p6_1, filename_n0p7_0, filename_n0p7_1]
ap_list = [ap_n0p6_0, ap_n0p6_1, ap_n0p7_0, ap_n0p7_1]

# For better plot visualization we are going to duplicate the bottom row and place it on top, and the leftmost column and place it on the right (periodic lattice)
Nx = 512 + 1
Ny = 512 + 1
Lx = 500.0
Ly = 500.0

x=np.arange(0.0, Nx)*(Lx/(Nx-1))
y=np.arange(0.0, Ny)*(Ly/(Ny-1))

xx=np.zeros((Nx,Ny)) 
yy=np.zeros((Nx,Ny)) 
for i in range(0,Nx): 
    for j in range(0,Ny): 
        xx[i,j]=x[i] 
        yy[i,j]=y[j] 

for i_ap in range(4):
    ap = ap_list[i_ap]
    ap_expanded = np.zeros((Nx,Ny))

    # Duplicate the bottom and left data
    for i in range(0,Nx):
        for j in range(0,Ny):
             ap_expanded[i,j] = ap[i % (Nx-1), j % (Ny-1)]
    filename = filename_list[i_ap]

# print a little information so we can check what is in the files
    print("i_ap, filename:", i_ap, " ", filename)
    print("x[0,0], y[0,0], ap[0,0]: ", xx[0,0], yy[0,0], ap[0,0])
    print("x[1,0], y[1,0], ap[1,0]: ", xx[1,0], yy[1,0], ap[1,0])
    print("x[0,1], y[0,1], ap[0,1]: ", xx[0,1], yy[0,1], ap[0,1])
    print(" ")

    fig_contour, fig_contour_axes = plt.subplots(nrows = 1, ncols = 1, figsize=(6, 5), dpi=300, tight_layout=True)
    cmap = plt.get_cmap('Reds')
    levels = np.linspace(np.min(ap), np.max(ap), 20)
    contour_levels = levels[::3]
    fig_contour_contourf = fig_contour_axes.contourf(xx, yy, ap_expanded, levels=levels, cmap=cmap)
    fig_contour_contour = fig_contour_axes.contour(xx, yy, ap_expanded, levels=contour_levels, colors='black', linewidths=0.5)
    cbar = fig_contour.colorbar(fig_contour_contourf, location='right', ticks=contour_levels)
    cbar.ax.yaxis.set_major_formatter(FormatStrFormatter('%.3f'))
    cbar.set_label("a: predator density", fontsize=14)
    fig_contour_axes.set_aspect('equal')
    fig_contour_axes.set_xlabel(r'$x/2\pi$', fontsize=14)
    fig_contour_axes.set_ylabel(r'$y/2\pi$', fontsize=14)

    fig_contour.tight_layout()
    fig_contour.savefig(filename + '.pdf', bbox_inches='tight', transparent=True) 
    plt.show()
    plt.close(fig_contour)

