Normalisation of a spectra and display of data (matlab to python)

Question

I have a spectra, the code was originally matlab and I have converted it to python. There are 10 observations in this example, 29 orders. It has been continuum normalised, however the plot looks very messy and I think it has something to do with the indexing.

enter image description here. this is the matlab plot what it should look like

enter image description here this is my python plot

This is the current code

import os
import numpy as np
import h5py
from scipy.io import loadmat
import matplotlib.pyplot as plt
%matplotlib inline
%matplotlib widget

# Directory where the MATLAB files are located
directory = '/Users/xxxxxxx/Desktop/Tau Ceti/'

# List all files in the directory matching the pattern '202*.mat'
file_list = [file for file in os.listdir(directory) if file.startswith('202') and file.endswith('.mat')]

# Initialize empty lists or arrays to store the data
ALLcalspec1 = []
ALLcalspec2 = []
ALLfib4spec = []
ALLfib4spec0 = []
ALLfib4specwave = []
ALLfib4specbNc = []
ALLfib4specerr = []
ALLfib4specS = []
ALLJD = []

# Loop through each .mat file and load the data
for s, mat_file in enumerate(file_list):
    data = loadmat(os.path.join(directory, mat_file))
    
    ALLcalspec1.append(data.get('calspec1', np.array([])))
    ALLcalspec2.append(data.get('calspec2', np.array([])))
    ALLfib4spec.append(data.get('fib4_spec', np.array([])))
    ALLfib4spec0.append(data.get('fib4_spec0', np.array([])))
    ALLfib4specwave.append(data.get('fib4_spec_wave', np.array([])))
    ALLfib4specbNc.append(data.get('fib4_spec_bNc', np.array([])))
    ALLfib4specerr.append(data.get('fib4_specerr', np.array([])))
    ALLfib4specS.append(data.get('fib4_specS', np.array([])))
    ALLJD.append(data.get('JD', np.array([])))

# Directory where the synthetic spectrum MATLAB file is located
synthetic_spectrum_file = '/Users/xxxxxxx/Desktop/Tau Ceti/synth_teff4500_logg45.mat'

# Load the synthetic spectrum data using scipy's loadmat function
synthetic_spectrum_data = loadmat(synthetic_spectrum_file)

# Extract the relevant data from the loaded dictionary
synthetic_wavelengths = synthetic_spectrum_data.get('wavelength', np.array([]))
synthetic_intensities = synthetic_spectrum_data.get('intensity', np.array([]))

# Load the 'bcorr_fib4_all.mat' file using h5py
bcorr_file_path = '/Users/xxxxxx/Desktop/Tau Ceti/bcorr_fib4_all.mat'

with h5py.File(bcorr_file_path, 'r') as bcorr_file:
    # List the available keys or group names in the file
    available_keys = list(bcorr_file.keys())

# Assuming you have already loaded ALLfib4spec as a NumPy array with shape (29, 10, num_files)
# Where 'num_files' is the number of data files you have

# Permute the dimensions to match the MATLAB code
permuted_ALLfib4spec = np.transpose(ALLfib4spec, (2, 1, 0))  # or np.swapaxes(ALLfib4spec, 0, 2)

# Extract the 15th order flux for all observations
order_15_flux = permuted_ALLfib4spec[14, :, :]  # 14 because Python uses 0-based indexing

# Calculate the median flux for order 15 across all observations
median_flux_order_15 = np.median(order_15_flux, axis=0)

# Define the threshold for a "bad" observation
threshold = 10000  # 10000 electrons in flux for order 15

# Create the 'bad' boolean mask
bad = median_flux_order_15 < threshold

# Create an index list of observations that are not bad
good_indices = [i for i, is_bad in enumerate(bad) if not is_bad]

# Filter the NumPy arrays based on the 'good_indices' list
ALLfib4spec = [ALLfib4spec[i] for i in good_indices]
ALLfib4spec0 = [ALLfib4spec0[i] for i in good_indices]
ALLfib4specbNc = [ALLfib4specbNc[i] for i in good_indices]
ALLfib4specerr = [ALLfib4specerr[i] for i in good_indices]
ALLfib4specwave = [ALLfib4specwave[i] for i in good_indices]
ALLJD = [ALLJD[i] for i in good_indices]
ALLcalspec1 = [ALLcalspec1[i] for i in good_indices]
ALLcalspec2 = [ALLcalspec2[i] for i in good_indices]

# Make a copy of ALLfib4specbNc to store continuum-corrected data
nALLfib4specbNc = ALLfib4specbNc.copy()

# Iterate through observations
for obs in range(len(nALLfib4specbNc)):
    print(f'obs: {obs + 1}', end=' ')
    
    # Iterate through orders
    for ord in range(len(nALLfib4specbNc[obs])):
        x = nALLfib4specbNc[obs][ord]
        tx = np.arange(len(x))
        
        # Fit a straight line (first-degree polynomial) to the data
        p = np.polyfit(tx, x, 1)
        
        # Calculate the continuum-corrected data
        xf = np.polyval(p, tx)
        nA = nALLfib4specbNc[obs][ord] / xf
        
        # Update nALLfib4specbNc with the continuum-corrected data
        nALLfib4specbNc[obs][ord] = nA
        
#Now nALLfib4specbNc contains the continuum-corrected data

# Convert nALLfib4specbNc to a NumPy array
nALLfib4specbNc = np.array(nALLfib4specbNc)

# Threshold values
upper_threshold = 1.2
lower_threshold = 0.01

# Apply upper and lower threshold to nALLfib4specbNc
nALLfib4specbNc = np.where(nALLfib4specbNc > upper_threshold, upper_threshold, nALLfib4specbNc)
nALLfib4specbNc = np.where(nALLfib4specbNc < lower_threshold, lower_threshold, nALLfib4specbNc)

# Save the data
prepared_data = {
    'nALLfib4specbNc': nALLfib4specbNc,
    'ALLcalspec1': ALLcalspec1,
    'ALLcalspec2': ALLcalspec2,
    'synthetic_wavelengths': synthetic_wavelengths,
    'synthetic_intensities': synthetic_intensities,
    'bcorr_fib4_all': bcorr_fib4_all
}

# Save the prepared data to a .mat file
savemat('Prepared.mat', prepared_data)

# Set the figure size (adjust these values as needed)
plt.figure(figsize=(14, 6))  # Adjust width and height as desired

# Iterate through observations
for obs in range(len(nALLfib4specbNc[0, 0, :])):  # Using the shape of nALLfib4specbNc to get the number of observations
    # Iterate through orders
    for ord in range(len(nALLfib4specbNc)):  # Using the shape of nALLfib4specbNc to get the number of orders
        
        # Extract the wavelength and spectrum data for the current observation and order
        wavelength_obs_ord = ALLfib4specwave[ord][:, obs]
        spectrum_obs_ord = nALLfib4specbNc[ord][:, obs]
        
        # Transpose the data to match MATLAB's behavior
        wavelength_obs_ord = wavelength_obs_ord.T
        spectrum_obs_ord = spectrum_obs_ord.T
        
        # Plot the order
        plt.plot(wavelength_obs_ord, spectrum_obs_ord, label=f'Order {ord + 1}')

# Add labels and title
plt.xlabel('Wavelength')
plt.ylabel('Spectrum')
plt.title('Spectrum after Normalization (All Observations and Orders)')
plt.grid(True)

# Show the plot
plt.show()