simulationfunctions.py

from astropy.io import fits
import matplotlib.pyplot as plt
import numpy as np
from astropy.nddata.utils import Cutout2D
from astropy.wcs import WCS
from astropy import table
import pandas as pd
import math
from math import pi
from PIL import Image
from numpy import asarray
import yaml
from congrid import congrid
import generate_PSF
import numpy as np
import scipy.interpolate
import scipy.ndimage
import glob
import matplotlib.animation as animation
import imageio
import imageio.v3 as iio
import os
from scipy import interpolate
from matplotlib.colors import LogNorm
from scipy import signal
import scipy.constants as sc
import pdb
import copy 
import skimage.io
from skimage.util import random_noise
from datetime import datetime


# ===========================================================================
# :: PyISH: Python Integral Field Spectrograph Simulation ::
# ===========================================================================
# Below is every function in this simulation! 




# [ Image Processing ] -----------------------------------------------
# Image processing reads in the excel files you filled out specifiying your paths to the spatial maps, spectrum, specs, etc. 
# It will seperate the spectrum into flux contributions from the continuum vs flux from each emission
# It will then use each flux contribution to scale either the emission spatial maps, or the continuum spatial maps
# It will also convolve each map with the respective PSF
# Overall, image processing creates the input 3D cubes, that are scaled to Rayleighs and convolved with the PSF
# -----------------------------------------------------------------------------
def imageprocessing(widthof1slice, bandwidth, detectorbandpass, bandpassstart, detectorFOV):
	
	#Reading excel file
	sim_inputs = pd.read_excel('inputdataIFUSim.xlsx')
	path_to_spectrum = sim_inputs.loc[0, 'path_to_spectrum']

	#Pulls continuum info, already scaled to rayleighs
	albedowavelength, albedoflux = np.loadtxt(path_to_spectrum, skiprows=1, usecols=(0, 1), delimiter=',', unpack = True)


	#reading in inputs from excel file 
	lines = sim_inputs['Ion'].tolist()
	emissionline_wl = sim_inputs['Wavelength (nm)'].tolist()
	widthemissionlines = sim_inputs['delta_lambda (nm)'].tolist()
	
	# mask with all true values 
	mask = np.full_like(albedowavelength, True, dtype=bool)


	#Fun little loop that masks the emission lines, them fits a line to the remaining data
	#this remaining fit represents the flux contribution from the continuum

	for peak, width in zip(emissionline_wl, widthemissionlines):

		mask &= ~(((peak-(width/2)) <= albedowavelength) & (albedowavelength < (peak+(width/2))))


	#Essentially if you have a spectrum like this

	#  Y-Axis
	#   ^
	#   |          *
	#   |        *   *
	#   |       *     *                  *
	#   |      *       *               *   *
	#   |     *         *             *     *
	#   |    *           *           *       *
	#   |  *               * *  *  *           *
	#   |*                                      *
	#   +-------------------------------------------> X-Axis 
	#                       

	#Essentially if you have a spectrum like this 

	emissionsubtracted_spectrum_wl = albedowavelength[mask]
	emissionsubtracted_spectrum_flux = albedoflux[mask]


	#The mask will do this:

	#  Y-Axis
	#   ^
	#   |           
	#   |             
	#   |                                
	#   |                                 
	#   |                               
	#   |    *                            
	#   |  *                * *  *  *           *
	#   |*                                       *
	#   +-------------------------------------------> X-Axis 
	#                       

	
	
	
	#gets a function of the data
	#if I change fill value to extrapolate, I get a warning (I think it just becomes 0)
	function_spectrum = interpolate.interp1d(albedowavelength, albedoflux, fill_value= np.nan, bounds_error=False)

	function_continuum = interpolate.interp1d(emissionsubtracted_spectrum_wl,emissionsubtracted_spectrum_flux, fill_value= 'extrapolate', bounds_error=False)
	
	


	#Then the intepolation fills in. This is what your "function_continuum" would look like in your example:

	#  Y-Axis
	#   ^
	#   |           
	#   |             
	#   |                                
	#   |                                 
	#   |                               
	#   |    *  -  -  -                      
	#   |  *             -   * *  *  *  -  - - *
	#   |*                                       *
	#   +-------------------------------------------> X-Axis 
	
	#and your "function_spectrum is just a fit of your original spectrum "






	#This loop will open the continuum spatial map, then go through each wavelength in youur bandpass and flux-calibrate the continuum map to the flux at that wavelength step
	#In addition, when the loop reaches a specified emission line, it will flux-calibrate the emission to the flux contribution from just that emission line (total flux - continuum flux = emission flux contribution)
	#Calibrates to the center of each bin (wavelength), assuming 1 bin per bandwidth

	#NOTE: If you don't have a high resolution map of your source for the continuum (like a surface reflectance map), for your "continuum map" user input, input an array of 0's matching the shape of your emission line maps
	#PyISH will still work

	contdata_listR = []
	continuummap = sim_inputs.loc[0, 'path to continuum map']
	with fits.open(continuummap) as hdul:
			contdata = hdul[0].data

			
			#sets values in image between 0-1
			contdata = np.abs(contdata/np.max(contdata))
			

			#k is what goes through bandpass
			k = bandpassstart + bandwidth
			index = 0

		
			
			emission_wavelength_dict = [[] for _ in range(0,len(lines))]
			emission_flux_dict = [[] for _ in range(0,len(lines))]

			
				
			
			while k <= bandpassstart + detectorbandpass:

				#print(k)

				#figuring out where the center of the wavelength bin is.
				#Example: IF there is 1 data point between 100nm and 101nm, lambda correspondence would be 100.5nm
				lambdacorrespondence = k - (bandwidth/2)
				

				
				spectrumflux = function_spectrum(lambdacorrespondence)
				

				if spectrumflux < 0 or math.isnan(spectrumflux): 
					spectrumflux = np.float64(0) 
					

				
				#Figures out flux contribution from the continuum at the specified wavelength 
				continuumflux = function_continuum(lambdacorrespondence)

				

				if continuumflux < 0 or math.isnan(continuumflux): 
					continuumflux = np.float64(0)
					


				#flux calibrates the continuum spatial map to the flux at that wavelength
				scaledcontdata = contdata * continuumflux

				

				#saves all of these bad boys in a list 
				contdata_listR.append(scaledcontdata)
				

				#loop inside of a loop (fancy) that checks if we have reached a emission line
				#if k is the wavelength of a listed emission line, then we flux calibrate the emission line map as well! 
				ion_counter = 0
				for ion, peak, width in zip(lines, emissionline_wl, widthemissionlines):
					#determines the range for an emission line
					lower = (peak-(width/2))
					upper = (peak+(width/2))


					if lower <= k < upper:
						#list of the wavelengths that are associated with emission lines
						emission_wavelength_dict[ion_counter].append(k)
						#if in the range, flux calibrate the associated map
						em_flux = spectrumflux - continuumflux
					
						if em_flux < 0:
							em_flux = 0 
						#list of the fluxes associated with emission lines
						emission_flux_dict[ion_counter].append(em_flux)

						
						
						break
					ion_counter+=1
						
						
			


				k += bandwidth
				index += 1
	
	#Just to recap on the most recent loop
	#the loop goes across the x-axis (wavelength) in steps consistent with the resolution
		#  Y-Axis (Amplitude)
	#   ^
	#   |          * 
	#   |        *   * 
	#   |       *     *                                  *
	#   |      *       *                               *   *
	#   |     *         *                             *     *
	#   |    *           *                           *       *
	#   |  *               *                       *           *
	#   |*                   * * * * * * * * * * *               *
	#   +----------|----------------|-----------------------------> X-Axis (wavelength)
	#   0          |                |           
	#              ^                ^
	#   [ when k gets here ]          [ when k gets here ]
	# flux calibrate emission         flux calibrate the continuum spatial map only!
	# total - continuum = emission
	# flux calibrate continuum
	# with continuum function




	#All of this in between the hashtags is how I made one of the figures in my paper
	########################################################################################
	#y1 = function_spectrum(albedowavelength)
	#y2 = function_continuum(albedowavelength)
	#plt.plot(albedowavelength, y1, color = "Blue", label = "UV Spectrum of Europa")	
	#plt.plot(albedowavelength, y2, color = "Red", label = "Continuum of Europa", linestyle='dashed')

	#plt.fill_between(albedowavelength, y2, color = "Purple", label = "Continuum of Europa", linestyle='dashdot')
	#plt.grid()
	#plt.xlabel("Wavelength [nm]")
	#plt.ylabel("Flux [R]")
	#plt.legend()
	#plt.yscale('log')
	#plt.savefig('spectrumeuropa.png', dpi = 600)
	#plt.show()	

	#commented out in case you wanna pause here
	#pdb.set_trace()
	#########################################################################################
	
	
	model_files = sim_inputs['path to model spectral map'].tolist()



	inputs = pd.read_excel('InputInformationIFU.xlsx')
	pixelscale = inputs.loc[0, 'Pixel_Scale'] #mas/pixel
	eac_files = inputs['path to eac yamls'].tolist()

	eac = inputs.loc[0, 'EAC']
	


	
	data_list = [[] for _ in range(0,len(model_files))]

	print("step size: ", bandwidth, " nm"  )
	psf_option = input("Option 1: Generate a new PSF at every wavelength step, or Option 2: Use one wavelength - specified PSF for the entire bandpass? Type `1` for Option 1, `2` for Option 2. ")

	if psf_option == '2':
		psf_option_wl = float(input("You have selected Option 2. What wavelength would you like the PSF? Please specify in nm: "))

		psf_option2 = generate_PSF.create_psf(psf_option_wl, eac_files[int(eac) - 1], pixelscale, FOV = np.array([detectorFOV, detectorFOV]), EAC = str(int(eac)), display = False)


	counter = 0
	#This next loop is applying the appropriate PSF to every flux - calibrated map we just did
	for model, ion_list, wavelengthlist, fluxlist in zip(model_files, lines, emission_wavelength_dict, emission_flux_dict):



		with fits.open(model) as hdul_1:
			modeldata = hdul_1[0].data

			#my images were flipping and rotating for some reason so I have these commented out just in case it happens to you
			#modeldata = np.rot90(modeldata, k=2)

			modeldata = modeldata / (np.max(modeldata))

			#modeldata = np.flipud(modeldata)

			#pdb.set_trace()
			



			
			
		

			for j in range(0, len(wavelengthlist)):
					
					
					if psf_option == '1':

						psf = generate_PSF.create_psf(wavelengthlist[j], eac_files[int(eac) - 1], pixelscale, FOV = np.array([detectorFOV, detectorFOV]), EAC = str(int(eac)), display = False)
						

					if psf_option == '2':
						psf = psf_option2
						
					
					#convolving PSF
					modelwpsf = signal.convolve2d(modeldata , psf, mode = 'same')

					#if you don't want to convolve with the psf at all, comment out line 350 (where I define modelwpsf, and switch "modelwpsf" with 'modeldata' in the line below. (you can see that in my commented version))
					#convolving takes a lot of time, so this helps check things 
					#scaledmodeldata = fluxlist[j] * modeldata
					scaledmodeldata = fluxlist[j] * modelwpsf
					
					#uncomment this if you want to see what your images look like post - convolving. beware, this will go through your entire bandpass in tiny steps
					#plt.imshow(scaledmodeldata)
					#plt.show()


						

					data_list[counter].append(scaledmodeldata)
			time = datetime.now()
			e = time.strftime("%H:%M:%S")	
			print("just got through image processing ion: " , ion_list," current time:",e)			



			counter += 1
			
	
						
	time1 = datetime.now()	
	f =  time.strftime("%H:%M:%S")	
	print("just got through image processing completely, current time:", f)	
	return data_list, albedowavelength, albedoflux, contdata_listR, emission_wavelength_dict, widthemissionlines, emissionline_wl, lines



# [Telescope Calibration] -----------------------------------------------
# This will take your input 3D cube that you made in the image processing function, and convert the Rayleighs to a photon count 
# Rayleighs -> photons depends on instrument specs like efficiencies, exposure time, etc
# which is why these are all user inputs!
# 
# -----------------------------------------------------------------------------

def telescopecalibration(data_list, emission_wavelength_dict, contdata_listR, bandpassstart, detectorbandpass, bandwidth, widthemissionlines, emissionline_wl, detectorpixelscale, splitdetector, cutoutsonxaxis, cutoutsonyaxis ):

	
	
	#extracting efficiencies 
	sim_inputs = pd.read_excel('inputdataIFUSim.xlsx')

	arcsecdata = sim_inputs.loc[0, 'arcsecspectra']

	


	pixelsasdetector = int(arcsecdata / detectorpixelscale)
	heightof1slice = splitdetector[cutoutsonyaxis-1][cutoutsonxaxis-1].shape[0]



	path_to_grating = sim_inputs.loc[0, 'path to grating yaml']

	#these are pulling all the efficiencies from the yaml files
	#I could have made this much more efficient but I am ~lazy~ 
	with open(path_to_grating, 'r') as file:
		grating_efficiency = yaml.safe_load(file)

	grating = np.array([grating_efficiency['wavelength'], grating_efficiency['grating_efficiency']])
	#2d numpy array out of yaml file
	x_g,y_g = grating

	path_to_QE = sim_inputs.loc[0, 'path to QE yaml']
	with open(path_to_QE, 'r') as file:
		QE = yaml.safe_load(file)

	MCP_QE = np.array([QE['wavelength'], QE['QE']])
	x_QE,y_QE = MCP_QE

	path_to_coating = sim_inputs.loc[0, 'path to coating yaml']
	with open(path_to_coating, 'r') as file:
		coating_efficiency = yaml.safe_load(file)

	coating = np.array([coating_efficiency['wavelength'], coating_efficiency['reflectivity']])
	x_c,y_c = coating


	inputs = pd.read_excel('InputInformationIFU.xlsx')

	
	reflections = inputs.loc[0, 'reflections'] 
	
	
	Inscribed_diameter = inputs.loc[0, 'Inscribed_diameter'] #m

	AGeo = pi * ((Inscribed_diameter/2)**2)

	#All of this commented out is how I made the efficiency plots for the .yaml files
	'''
	
	plt.plot(x_g,y_g, color = "Blue", label = "Grating Efficiency")

	plt.plot(x_QE,y_QE, color = "Green", label = "Quantum Efficiency",linestyle='dashed')

	plt.plot(x_c,y_c, color = "Purple", label = "Coating Reflectivity", linestyle='dashdot')

	plt.xlim(90,400)
	plt.ylim(0,1)
	plt.xlabel('Wavelength [nm]')
	plt.legend()
	plt.ylabel('Efficiency')
	#plt.grid()
	#plt.show()
	plt.savefig('efficiencies.png', dpi = 600)
	#plt.show()
	pdb.set_trace()

	'''

	#Exposure Time
	ExpTime = inputs.loc[0, 'ExpTime']

	#pixel scale (converted to radians)
	pixelscale_r = inputs.loc[0, 'Pixel_Scale'] * 0.001 * (pi/(180*3600))

	domega  = pixelscale_r * pixelscale_r
	


	


	#continuum scaling to photon count 	
	contdata_listphoton = []
	scaledHRIdata = [[] for _ in range(0,len(data_list))]
	indexingcontinuum = 0 
	k = bandpassstart + bandwidth
	
	ep = []
	test = []
	while k <= (bandpassstart + detectorbandpass):
			
		#sets scaling to middle of the bin
		lambdacorrespondence = k - (bandwidth/2)
		
				
		
		difference_g_cont = np.absolute(x_g - k)
		indicies_g_cont = difference_g_cont.argmin() 
		grating_efficiency_continuum = float(y_g[indicies_g_cont])
			

		difference_QE_cont = np.absolute(x_QE - k)
		indicies_QE_cont = difference_QE_cont.argmin() 
		detector_loss_continuum = float(y_QE[indicies_QE_cont])
			
			

		difference_c_cont = np.absolute(x_c - k)
		indicies_c_cont = difference_c_cont.argmin() 
		coating_efficiency_cont = float(y_c[indicies_c_cont])
			

		epsilon_cont = grating_efficiency_continuum*detector_loss_continuum*(coating_efficiency_cont**reflections)

		ep.append(epsilon_cont)



		AEff_lamb = AGeo * epsilon_cont


		#this below is the equation for Rayleighs -> photons. It depends on effective area (which is geometric area * efficiencies), exposure time, pixel scale squared, and some constants

		contdata_photoncount = contdata_listR[indexingcontinuum] *  AEff_lamb * ExpTime * domega *((1e10)/(4*pi))


		#the couple of lines below are reshaping your maps!
		#the input spatial maps are not neccesarily representative of how the soruce would be captured by an HWO IFS
		#more specifically, the pixel pitch needs adjustment. If your IFS has a field of view of 5 arcseconds, and a pixel pitch of 0.1 arcseconds/pixel,
		#these next couple lines makes sure that each image is 50 x 50 pixels 


		#first reshapes to the correct amount of pixels tall 
		contdata_photoncount = resample(contdata_photoncount, (pixelsasdetector, contdata_photoncount.shape[1]), method = 'linear')

		#then adds extra 0 space to match detector height
		#if your image was too little pixels, extra pixels will be filled in as 0's (as if no data was captured)
		#if your image was too many pixels, iteratively chops off top and bottom of image, one pixel at a time to reach desired size
		contdata_photoncount = rebin(contdata_photoncount,  heightof1slice)




		#Next two lines do the exact same thing on the x-axis
		contdata_photoncount = resample(contdata_photoncount, (contdata_photoncount.shape[0], pixelsasdetector), method = 'linear')

		contdata_photoncount = rebinx(contdata_photoncount, heightof1slice)
		






		contdata_listphoton.append(contdata_photoncount)

		

		#the loop above went from rayleighs -> photons, and made sure the sizes were coreect for the continuum spatial map images
		#the loop below does the exact same thing for the emission lines
		#again I could make this more efficient but as I always say...if it aint broke dont fix it
		emissioncounter = 0
		
		for peak, width, wavelengthlist in zip(emissionline_wl, widthemissionlines, emission_wavelength_dict):

					lower = (peak-(width/2))
					upper = (peak+(width/2))


					if lower <= k < upper:

						
						values = np.absolute(wavelengthlist - k)
						indexofinterest = np.argmin(values)
						
						#scales to photon count
						image_photoncount = data_list[emissioncounter][indexofinterest] *  AEff_lamb * ExpTime * domega *((1e10)/(4*pi))

						

						#first reshapes to pixels as detector high
						image_photoncount = resample(image_photoncount, (pixelsasdetector, image_photoncount.shape[1]), method = 'linear')

						#then adds extra 0 space to match detector height
						image_photoncount = rebin(image_photoncount,  heightof1slice)


						#first reshapes to pixels as detector wide
						image_photoncount = resample(image_photoncount, (image_photoncount.shape[0], pixelsasdetector), method = 'linear')

						#then matches detector width
						#height of one slice is in pixels, this sets fov in input transform
						#make sure this is actually true 
						image_photoncount = rebinx(image_photoncount, heightof1slice)

						
						
						
						

						scaledHRIdata[emissioncounter].append(image_photoncount)

						
						
						

						


						
						break

					emissioncounter += 1 


		

		k += bandwidth
		indexingcontinuum = indexingcontinuum + 1

	
	
	
	time1 = datetime.now()	
	e = time1.strftime("%H:%M:%S")	
	print("just got through telescope calibration completely, current time:", e)	
	return  scaledHRIdata, contdata_listphoton

 

# [Image resampling methods] -----------------------------------------------
#The functions below are called within the modules of PyISH, not independently

# -----------------------------------------------------------------------------


__all__ = ['resample', 'reshape_image_to_4d_superpixel']



def resample(orig, dimensions, method='linear', center=False, minusone=False):
    """
    Returns a new `numpy.ndarray` that has been resampled up or down.

    Arbitrary resampling of source array to new dimension sizes.
    Currently only supports mai ntaining the same number of dimensions.
    To use 1-D arrays, first promote them to shape (x,1).

    Uses the same parameters and creates the same coordinate lookup points
    as IDL's ``congrid`` routine (which apparently originally came from a
    VAX/VMS routine of the same name.)

    Parameters
    ----------
    orig : `numpy.ndarray`
        Original input array.
    dimensions : `tuple`
        Dimensions that new `numpy.ndarray` should have.
    method : {``"nearest"``, ``"linear"``, ``"spline"``}, optional
        Method to use for resampling interpolation.
            * nearest and linear - Uses "n x 1D" interpolations calculated by
              `scipy.interpolate.interp1d`.
            * spline - Uses `scipy.ndimage.map_coordinates`
    center : `bool`, optional
        If `False` (default) the interpolation points are at the front edge of the bin.
        If `True`, interpolation points are at the centers of the bins
    minusone : `bool`, optional
        For ``orig.shape = (i,j)`` & new dimensions ``= (x,y)``, if set to `False`
        (default) ``orig`` is resampled by factors of ``(i/x) * (j/y)``,
        otherwise ``orig`` is resampled by ``(i-1)/(x-1) * (j-1)/(y-1)``.
        This prevents extrapolation one element beyond bounds of input array.

    Returns
    -------
    out : `numpy.ndarray`
        A new `numpy.ndarray` which has been resampled to the desired dimensions.

    References
    ----------
    https://gcc02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fscipy-cookbook.readthedocs.io%2Fitems%2FRebinning.html&data=05%7C02%7Cgrace.m.sweetak%40nasa.gov%7C95cf4e9028f8472607d808ddb807fbfe%7C7005d45845be48ae8140d43da96dd17b%7C0%7C0%7C638869066600300365%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=LZeINPbSOmvBclh8v5BW8VSpLb%2BtzEMOplAWlvzgYVo%3D&reserved=0
    """

    # Verify that number dimensions requested matches original shape
    if len(dimensions) != orig.ndim:
        raise UnequalNumDimensions("Number of dimensions must remain the same "
                                   "when calling resample.")

    # TODO: Will this be okay for integer (e.g. JPEG 2000) data?
    if orig.dtype not in [np.float64, np.float32]:
        orig = orig.astype(np.float64)

    dimensions = np.asarray(dimensions, dtype=np.float64)
    m1 = np.array(minusone, dtype=np.int64)  # array(0) or array(1)
    offset = np.float64(center * 0.5)       # float64(0.) or float64(0.5)

    # Resample data
    if method in ['nearest', 'linear']:
        data = _resample_nearest_linear(orig, dimensions, method,
                                        offset, m1)
    elif method == 'spline':
        data = _resample_spline(orig, dimensions, offset, m1)
    else:
        raise UnrecognizedInterpolationMethod(f"Unrecognized interpolation method requested: {method}")

    return data




def _resample_nearest_linear(orig, dimensions, method, offset, m1):
    """
    Resample Map using either linear or nearest interpolation.

    Parameters
    ----------
    orig : array-like
        Original data.
    dimensions : `tuple`
        Dimensions of resampled data.
    method : `str`
        Interpolation method passed to `~scipy.interpolate.interpn`
    offset : `float`
        Either 0 or 0.5, depending on whether interpolation is at the edge or
        centers of bins.
    m1 : 0 or 1
        For ``orig.shape = (i,j)`` & new dimensions ``= (x,y)``, if set to `False`
        (default) ``orig`` is resampled by factors of ``(i/x) * (j/y)``,
        otherwise ``orig`` is resampled by ``(i-1)/(x-1) * (j-1)/(y-1)``.
        This prevents extrapolation one element beyond bounds of input array.
    """
    old_coords = [np.arange(i, dtype=float) + offset for i in orig.shape]
    scale = (orig.shape - m1) / (dimensions - m1)
    new_coords = [(np.arange(dimensions[i], dtype=float) + offset) * scale[i] for i in
                  range(len(dimensions))]
    new_coords = np.stack(np.meshgrid(*new_coords, indexing='ij'), axis=-1)
    # fill_value = None extrapolates outside the domain
    new_data = scipy.interpolate.interpn(old_coords, orig, new_coords,
                                         method=method, bounds_error=False,
                                         fill_value=None)

    return new_data


def _resample_spline(orig, dimensions, offset, m1):
    """
    Resample Map using spline-based interpolation.
    """
    nslices = [slice(0, j) for j in list(dimensions)]
    newcoords = np.mgrid[nslices]

    newcoords_dims = list(range(newcoords.ndim))

    # make first index last
    newcoords_dims.append(newcoords_dims.pop(0))
    newcoords_tr = newcoords.transpose(newcoords_dims)

    # makes a view that affects newcoords
    newcoords_tr += offset

    deltas = (np.asarray(orig.shape) - m1) / (dimensions - m1)
    newcoords_tr *= deltas

    newcoords_tr -= offset

    return scipy.ndimage.map_coordinates(orig, newcoords)




def reshape_image_to_4d_superpixel(img, dimensions, offset):
    """
    Re-shape the two-dimensional input image into a four-dimensional array
    whose first and third dimensions express the number of original pixels in
    the "x" and "y" directions that form one superpixel. The reshaping makes it
    very easy to perform operations on superpixels.

    An application of this reshaping is the following. Let's say you have an
    array::

         x = np.array([[0, 0, 0, 1, 1, 1],
                       [0, 0, 1, 1, 0, 0],
                       [1, 1, 0, 0, 1, 1],
                       [0, 0, 0, 0, 1, 1],
                       [1, 0, 1, 0, 1, 1],
                       [0, 0, 1, 0, 0, 0]])

    and you want to sum over 2x2 non-overlapping sub-arrays. For example, you
    could have a noisy image and you want to increase the signal-to-noise ratio.
    Summing over all the non-overlapping 2x2 sub-arrays will create a
    superpixel array of the original data. Every pixel in the superpixel array
    is the sum of the values in a 2x2 sub-array of the original array.

    This summing can be done by reshaping the array::

         y = x.reshape(3,2,3,2)

    and then summing over the 1st and third directions::

         y2 = y.sum(axis=3).sum(axis=1)

    which gives the expected array::

         array([[0, 3, 2],
               [2, 0, 4],
               [1, 2, 2]])

    Parameters
    ----------
    img : `numpy.ndarray`
        A two-dimensional `numpy.ndarray` of the form ``(y, x)``.
    dimensions : array-like
        A two element array-like object containing integers that describe the
        superpixel summation in the ``(y, x)`` directions.
    offset : array-like
        A two element array-like object containing integers that describe
        where in the input image the array reshaping begins in the ``(y, x)``
        directions.

    Returns
    -------
    A four dimensional `numpy.ndarray` that can be used to easily create
    two-dimensional arrays of superpixels of the input image.
    """
    # make sure the input dimensions are integers
    dimensions = [int(dim) for dim in dimensions]

    # New dimensions of the final image
    na = int(np.floor((img.shape[0] - offset[0]) / dimensions[0]))
    nb = int(np.floor((img.shape[1] - offset[1]) / dimensions[1]))

    # Reshape up to a higher dimensional array which is useful for higher
    # level operations
    return (img[int(offset[0]):int(offset[0] + na * dimensions[0]),
                int(offset[1]):int(offset[1] + nb * dimensions[1])]).reshape(na, dimensions[0], nb, dimensions[1])




class UnrecognizedInterpolationMethod(ValueError):
    """
    Unrecognized interpolation method specified.
    """


class UnequalNumDimensions(ValueError):
    """
    Number of dimensions does not match input array.
    """








# [Continuum Slicer] -----------------------------------------------
#This takes the input 3D sube for just the continuum spatial maps (which are scaled to a photon count and correctly sizes according to the pixel scale)
# It cuts the images up into the amount of slices as defined by the user inputs
#and saves these "contcutouts" = continuum cutouts 
# -----------------------------------------------------------------------------

def continuumslicer(contdata_listphoton, numberIFUslices):
	

	contcutouts = [[] for _ in range(0,len(contdata_listphoton))]
	
	for j in range(0, len(contdata_listphoton)):

		
		#retrieve image size
		xlength = contdata_listphoton[j].shape[1]
		ylength = contdata_listphoton[j].shape[0]

		#xcenter and ycenter are the centers of each slice
		#its a little weird but it works
		#  +-------+-------+
		#  |       |       |
		#  |       |       |
		#  |       |       |
		#  |       |       |
		#  |       |       |
		#  |       |       |
		#  |       |       |
		#  +---|----+-------+
		#      |
		#      ^ 
		#if my image is 8 pixels across, and I want 2 IFS slices, then xcenter = 2
		#2 pixels to the center of each slice

		ycenter = ylength/(numberIFUslices*2)
		xcenter = xlength/(numberIFUslices*2)

			
		#each image needs to be split into 100 IFS slices, then j iterates to the next image
		for i in range(0,numberIFUslices):

			#cut the slice
			contcutout = Cutout2D(contdata_listphoton[j], (xcenter, ylength/2), (ylength, xlength/numberIFUslices))


					
			#save the data, move to the next center and repeat
			data_to_save = contcutout.data
			ycenter = ycenter + (ylength/numberIFUslices)
			xcenter = xcenter + (xlength/numberIFUslices)

		
			
			contcutouts[j].append(
				 data_to_save
				  )
			
			#this is how I made a movie showin the cuts across each slice for image number 400 (replace 400 if have a really small bandpass of something)
			'''
			if j ==400:
				plt.imshow(contdata_listphoton[j], origin = 'lower', cmap = 'Blues_r')
				contcutout.plot_on_original(color = 'gray')
				plt.savefig(f"continuumcutouts/{i}.png")
				plt.close()
			'''
			#cutouts[0][2] represents the 1st image (H1), and the 2nd cutout/slice
	time1 = datetime.now()
	e = time1.strftime("%H:%M:%S")	
	print("just got through continuum cutouts completely, current time:", e)	
	
	return contcutouts


		
		

# [Image resampling methods] -----------------------------------------------
#This is basically the same function as continuumslicer but for the emission lines
#it's a little more complicated becuase we have multiple images for multiple emsisions, so have to have an extra loop in there
#I am sure there is a way to combine these functions but I do not feel like it!
#But yeah, this will cut up each photon-calibrated spatial map for each emission line, according to the number of IFS sluces
#will save everything in a list
# ---------------------------------------------------------------------------

def ifuslicer(scaledHRIdata, numberIFUslices, lines):
	
	cutouts_dict = {}




	#start by making empty lists of things
	for emissionlinenames, HRIlistforeachline in zip(lines, scaledHRIdata):
	    cutouts_dict[emissionlinenames] = [[] for _ in range(len(HRIlistforeachline))]

	rebinned_cutouts_dict = copy.deepcopy(cutouts_dict)


	
	

	counter = 0

	#top loop goes through each emission line

	for j, ionname in zip(range(0, len(scaledHRIdata)), lines):

		#this loop goes through all of the spatial maps associated with that ion
		for k in range(0, len(scaledHRIdata[counter])):

		
		
		
			xlength = scaledHRIdata[j][k].shape[1]
			ylength = scaledHRIdata[j][k].shape[0]

			#These are the centers of each slice
			ycenter = ylength/(numberIFUslices*2)
			xcenter = xlength/(numberIFUslices*2)

				
			#cutting them out as last function did
			for i in range(0,numberIFUslices):

			
				cutout = Cutout2D(scaledHRIdata[j][k], (xcenter, ylength/2), (ylength, xlength/numberIFUslices))
				

				data_to_save = cutout.data
				ycenter = ycenter + (ylength/numberIFUslices)
				xcenter = xcenter + (xlength/numberIFUslices)
				
				#this will show you one IFS slice for every spatial map for one emission line
				#if you have less than 35 IFS slices, then change the i 
				'''
				if i == 35:
					plt.imshow(scaledHRIdata[j][k], origin = 'lower', cmap = 'Blues_r')
					cutout.plot_on_original(color = 'gray')
					plt.xticks([])
					plt.yticks([])
					plt.xlim(0, scaledHRIdata[j][k].shape[1] - 1)
					plt.ylim(0, scaledHRIdata[j][k].shape[0] - 1)
					plt.colorbar()
					plt.show()
					
				'''
				

				cutouts_dict[ionname][k].append(data_to_save
						  )

				


		counter +=1 
	
	
			
		#cutouts_dict["Lya"[0][2] represents lyman alpha emission band, 0th wavelength in that band (around 119nm), 3rd cutout from that wavelength

		
		
	
	time1 = datetime.now()
	e = time1.strftime("%H:%M:%S")		
	print("just got through slicing ions completely, current time:", e)	
	return cutouts_dict, rebinned_cutouts_dict




# [Rebin] -----------------------------------------------
#This function is called within the other modules 
#it will iteratively add 0's to the top and bottom of your 2ds array, or delete rows iteratively depending on your new shape preference
# ---------------------------------------------------------------------------

#Making my own IDL rebin function
#adds extra 0's or deletes rows of an array to reshape
def rebin(array, new_shape_y): 
	
	'''
	Rebin function, that takes a 2D array of any size, and makes it into a (new_shape_y, 1) sized array

	The function starts by summing up all the rows, so you have an array of size (original, 1)

	When the old size is larger than the new shape, it will cut off top and bottom values

	For example : (5,1) going to (3,1)

 
    +----+
   	| 1  |
    +----+
    +----+
   	| 3  |
    +----+
    +----+
   	| 5  |
    +----+
    +----+
   	| 7  |
    +----+
    +----+
   	| 9  |
    +----+


	will produce :

    
    +----+
   	| 3  |
    +----+
    +----+
   	| 5  |
    +----+
    +----+
   	| 7  |
    +----+
	
	Oppositely, going from a smaller array to a larger array will repeat the end values
	For example: (3,1) to (5,1)


	+----+
   	| 3  |
    +----+
    +----+
   	| 5  |
    +----+
    +----+
   	| 7  |
    +----+


    will produce:

    +----+
   	| 3  |
    +----+
    +----+
   	| 3  |
    +----+
    +----+
   	| 5  |
    +----+
    +----+
   	| 7  |
    +----+
    +----+
   	| 7  |
    +----+

   '''
	


	oldarray_x = array.shape[1]
	oldarray_y = array.shape[0]
	#xarray = np.sum(array, axis = 1, keepdims = True)
	#print(xarray)


	toggle = 1
	while oldarray_y > new_shape_y:
		if toggle % 2 == 0:
			array = np.delete(array, 0, axis = 0)
			#print(array)
			
		else:
			array = np.delete(array, oldarray_y - 1 , axis = 0)
			#print(array)
			
		toggle = toggle + 1
		oldarray_y = oldarray_y - 1 

	while oldarray_y < new_shape_y:
		if toggle % 2 == 0:
	
			array = np.insert(array, 0, np.zeros(len(array[1].shape)), axis = 0)
			#print(array)
		else:
			array = np.insert(array, oldarray_y , np.zeros(len(array[1].shape)), axis = 0)
			#print(array)

		toggle = toggle + 1
		oldarray_y = oldarray_y + 1



	return array


# [Rebinx] -----------------------------------------------
#This function is the same as the rebin function but for the x axis
#so iteratively adds columns full of 0's to the beginning/end of arrays, or takes them aeay 
# ---------------------------------------------------------------------------

def rebinx(array, new_shape_x): 
	
	


	oldarray_x = array.shape[1]
	oldarray_y = array.shape[0]
	#xarray = np.sum(array, axis = 1, keepdims = True)
	#print(xarray)


	toggle = 1
	while oldarray_x > new_shape_x:
		if toggle % 2 == 0:
			array = np.delete(array, 0, axis = 1)
			#print(array)
			
		else:
			array = np.delete(array, oldarray_x - 1 , axis = 1)
			#print(array)
			
		toggle = toggle + 1
		oldarray_x = oldarray_x - 1 

	while oldarray_x < new_shape_x:
		if toggle % 2 == 0:
	
			array = np.insert(array, 0, np.zeros(len(array[0].shape)), axis = 1)
			#print(array)
		else:
			array = np.insert(array, oldarray_x , np.zeros(len(array[0].shape)), axis = 1)
			#print(array)

		toggle = toggle + 1
		oldarray_x = oldarray_x + 1



	return array


# [Input Transform] -----------------------------------------------
#This function is one of the two WORKHORSES of PyISH
#This function handles all of the user input trades for everything that you input
#This specific function is described in Sweetak et. al 2026 (JATIS)
#There is an updated version described in an SPIE proceeding
#the website that goes along with that update can be found at hwopyish.streamlit.app

#
#outputs: tells you whether your inputs are physical given the detector space that you specified. If so, it will give you updated specs based on user - input preferences 
# ---------------------------------------------------------------------------

def inputtransform():

	
	#reading inputs
	detectorinputs = pd.read_excel('InputInformationIFU.xlsx')

	detectorpixelscale = detectorinputs.loc[0, 'Pixel_Scale'] * 1e-3 #arcsec/pixel

	minFOV = detectorinputs.loc[0, 'FOV'] #FOV - arcseconds

	R = detectorinputs.loc[0, 'R']

	bandpass = detectorinputs.loc[0, 'bandpass'] #nm

	bandpassstart = detectorinputs.loc[0, 'bandpass_start'] #nm

	numberIFUslices = int(detectorinputs.loc[0, 'numberIFUslices'])

	lamb = detectorinputs.loc[0, 'lambda'] #nm
	bandwidth = lamb/R
	
	


	ylengthdetector = detectorinputs.loc[0, 'detector_size2']
	xlengthdetector = detectorinputs.loc[0, 'detector_size1']

	detectorsize = np.empty([ylengthdetector, xlengthdetector])

	pixelpitch = detectorinputs.loc[0, 'pixelpitch']  /1000 #in mm


	#reads in detector size in mm, then converts to pixels 
	detectorsizeinpixelsx = int(xlengthdetector/ pixelpitch)
	detectorsizeinpixelsy = int(ylengthdetector/ pixelpitch)
	detectorsizeinpixels = np.empty((detectorsizeinpixelsy,detectorsizeinpixelsx))

	
	
	#starts by making even cuts in x and y 
	cutoutsonyaxis = int(np.sqrt(numberIFUslices))
	cutoutsonxaxis = int(numberIFUslices/cutoutsonyaxis) 

	#ensures we have the correct amount of boxes for IFS slices on detector
	#each IFS slice corresponds to one "box", or one section of the detector
	#this is because every IFS slice has distinct spectra that is placed in distinct regions on the detector. so we need to make those sections
	if (cutoutsonyaxis*cutoutsonxaxis) < numberIFUslices: 
		cutoutsonxaxis = cutoutsonxaxis + 1

	
	#Sets minimum width and height needed for each slice based on the user inputs
	#the x axis, or width for each section corresponds to the bandpass, so need enough room for that
	#the factor of 2 is for nyquist sampline
	#the height corresponds to the FOV for eac hslice
	minwidthof1slice =  np.ceil((bandpass/bandwidth) * 2)
	pixelheightperslice = np.ceil(minFOV/detectorpixelscale)

	#decides which dimension input by the user will serve as x vs. y 
	if (minwidthof1slice > pixelheightperslice) and (detectorsizeinpixelsx < detectorsizeinpixelsy):

		ylengthdetector = detectorinputs.loc[0, 'detector_size1']
		xlengthdetector = detectorinputs.loc[0, 'detector_size2']

		detectorsize = np.empty([ylengthdetector, xlengthdetector])
		detectorsizeinpixelsx = int(xlengthdetector/ pixelpitch)
		detectorsizeinpixelsy = int(ylengthdetector/ pixelpitch)
		detectorsizeinpixels = np.empty((detectorsizeinpixelsy,detectorsizeinpixelsx))
		print(" ")
		print(" ")
		print("the width of 1 slice takes up more space than the height of one slice, so rotating detector to maximize space")
		print(" ")
		print(" ")
		print("---------------")
		print("-             -                                            ----------------------  ")
		print("-             -                                            -                    -  ")
		print("-             -               ------>                      -                    -  ")
		print("-             -                                            -                    -  ")
		print("-             -                                            ----------------------  ")
		print("-             -")
		print("---------------")
		print("narrow and tall               ------>                         short and wide")
		print(" ")
		print(" ")



	#initializes the minimum detector size given even cuts on x and y
	# these next two lines answer the question: if we have a certain orientation of sections, and each section has a certain min size, whats the total minimum detector space we need?
	minxdetectorsize  = minwidthof1slice*cutoutsonxaxis #pixels
	minydetectorsize = pixelheightperslice * cutoutsonyaxis #pixels


	
	#if both constraints for x and y are too small, then we can't do any tradeoff with IFS slice orientation
	if minxdetectorsize > detectorsizeinpixels.shape[0] and minydetectorsize > detectorsizeinpixels.shape[1]:

			print("Your detector is too small given the spectral and spatial constraints.", "The detector size given was ", detectorsize.shape, 
				" mm, but to acheieve the bandpass, resolution, and FOV, you need a minimum detector of size", minydetectorsize*pixelpitch, ", ", minxdetectorsize*pixelpitch, "mm")
			#possible detector is variable that is passed to decide if the main function will continue. If it is not possible to have an IFS with your specific constraints, possible detector = 0
			possibledetector = 0

	

	#if minimum width is too large to fit on detector 
	#option B (which is this loop) is described in Sweetak et al 2026
	elif minxdetectorsize > detectorsizeinpixels.shape[1]:
		
		
		print('Right now, the organization of sections in the detector does not allow you to achieve the desired spectral resolution and bandpass. Attempting to adjust spacing on detector to fit slices with Option B...' )


		while minxdetectorsize > detectorsizeinpixels.shape[1]:


			#adjust shaping
			#this nested loop attempts to adjust how each section is oriented. 
			#it will iteratively remove a column, and add a row, making the sections wider but less tall

			#total amount of columns
			cutoutsonxaxis = cutoutsonxaxis - 1
			
			#total amount of rows
			cutoutsonyaxis = cutoutsonyaxis + 1
			

			#this check makes sure the total amount of sections still >= numberIFUslices
			IFUcheck = cutoutsonyaxis*cutoutsonxaxis


			while (IFUcheck) < numberIFUslices: 
		
				cutoutsonyaxis = cutoutsonyaxis + 1
				
				IFUcheck = cutoutsonyaxis*cutoutsonxaxis
				
				#if you can't remove any more columns, no adjustment can be made. we need to thus use a larger detector
				if cutoutsonxaxis == 0:
					print("Changing the orientation of the cutouts to achieve spectral constraints voided the spatial constraints needed.",
					" The only option to achieve both is to use bigger detector.")
					possibledetector = 0
					break
				

			
			#recalculate detector size needed given these cuts in the x and y axis 
			#minwidthof1slice and pixelheightperslice will never change - these are based on your bandpass, resolution, FOV, and pixel scale
			#but the number of cuts, or orientation of slices (rows vs columns) is changing, which will change your detector size
			#in this case, we are trying to minimumize the minxdetectorsize. That value will get smaller as minydetectorsize increased
			minxdetectorsize  = minwidthof1slice*cutoutsonxaxis
			minydetectorsize = pixelheightperslice * cutoutsonyaxis
			
			

			#This value "cutoutsonxaxis" may look misleading, because "0" could mean that there are 0 "cuts" to the detector, and we just have 1 column
			#BUT I use these variables in np.array_split, and so if you pass "1", that means just have 1 array (no splitting)
			#passing "0" results in a value error
			#so think of "cutoutsonxaxis" and "cutoutsonyaxis" as the total amount of columns and rows, respectively
			if cutoutsonxaxis == 0:
					
					possibledetector = 0
					break
			

			#in this if startement, the minydetectorsize got too large for what the user input is. 
			if minydetectorsize > detectorsizeinpixels.shape[0]:
					
				print("Changing the orientation of the cutouts to achieve spectral constraints voided the spatial constraints needed.",
				" The only option to achieve both is to use bigger detector.")
				possibledetector = 0


				break 
			
			#adjustment successful! possibledetector = 1, which means we can accomadate user inputs and the simulation continues
			if minydetectorsize <= detectorsizeinpixels.shape[0] and minxdetectorsize <= detectorsizeinpixels.shape[1] and cutoutsonxaxis != 0:
					print("After adjusting spacing, this detector size works for your given constraints! ")
					possibledetector = 1
					




	#if minimum height is too large to fit on detector 
	#this is option A, described in Sweetak et al 2026
	elif minydetectorsize > detectorsizeinpixels.shape[0]:



		print('Right now, the organization of sections in the detector does not allow you to achieve the desired FOV given the pixel plate scale. Attempting to adjust spacing on detector to fit slices with Option A...' )

		while minydetectorsize > detectorsizeinpixels.shape[0]:


			#adjust shaping
			#now, instead of removing a column and adding a row, we are removing a row, and adding a column
			#this makes sections taller and skinnier
			cutoutsonyaxis = cutoutsonyaxis - 1
			
			cutoutsonxaxis = cutoutsonxaxis + 1
			
			#checking that we still have enough sections
			while (cutoutsonyaxis*cutoutsonxaxis) < numberIFUslices: 
		
				cutoutsonxaxis = cutoutsonxaxis + 1
				
				#if we can make no more trades because there are no more rows, and still not enough sections, then cannot make the user inputs work
				if cutoutsonyaxis == 0:
					print("Changing the orientation of the cutouts to achieve spatial constraints voided the spectral constraints needed.",
					" The only option to achieve both is to use bigger detector.")
					possibledetector = 0
					break



			#recalculate detector size needed given these cuts in the x and y axis 
			#again, the minwidthof1slice and pixelheightperslice will never change, but the minxdetectorsize and minydetectorsize will 
			minxdetectorsize  = minwidthof1slice*cutoutsonxaxis
			
			minydetectorsize = pixelheightperslice * cutoutsonyaxis
			
			if cutoutsonyaxis == 0:
				
				possibledetector = 0
				break 
			
			#trading space on the y axis by adding space to the x axis (more columns) made the y-axis too big
			if minxdetectorsize > detectorsizeinpixels.shape[1]:
				print("Changing the orientation of the cutouts to achieve the desired FOV voided the spectral constraints needed.",
					" The only option to achieve both is to use bigger detector.")
				possibledetector = 0
				break 
			

			#the trades we made worked! 
			if minydetectorsize <= detectorsizeinpixels.shape[0] and minxdetectorsize <= detectorsizeinpixels.shape[1] and cutoutsonxaxis != 0:
				print("After adjusting spacing, this detector size works for your given constraints! ")
				possibledetector =1





	#if both height and width fit into the constraints
	#this is if no treades need to be made and the detector sizes are big enough with # rows =  # columns
	elif minydetectorsize <= detectorsizeinpixels.shape[0] and minxdetectorsize <= detectorsizeinpixels.shape[1]:
		print("Your detector is large enough for the spectral and spatial constraints!")
		possibledetector = 1

		
	#The detector is either possible or not possible. 
	#if not possible, simulation will not continue 
	#if detector setup is possible, it will split 2 empty arrays that are the same size as the detector
	#One empty array will be for the emission lines
	#the other empty array will be for the continuum
	#these will be added together in the end 


	if possibledetector == 0:
		print("Given the user inputs, I can not section the detector appropriately to fit all of the desired spectra.")
		print("So, I cannot continue. Either change the science constraints or provide a larger detector and rerun.")
		print("You are asking too much of me :(")
		print(" ")
		print(" ")
		print(" ")
		print(" ")
		print(" ")
		print(" ")
		print(" ")




	else:

		optimization = input("Would you like use the entire space avaliable on the given input detector? This would potentially widen your FOV and bandpass. If so, type y. Otherwise, type n to use the minimum detector size for your science constraints.")
		if optimization == "y":
			#assuming pixel format
			
			emptydetector = detectorsizeinpixels
		elif optimization == "n":
			#assuming pixel format
			#if you choose not to use your original input, it will use the minimum detector sizes determined above
			#these values are determined by the cutoutssonxaxis and cutoutssonyaxis, so if you have a square vs rectangular detector input,
			#these minimum detector sizes below may change depending on the orientation trades the simulation made
			minydetectorsize = int(minydetectorsize )
			minxdetectorsize = int(minxdetectorsize )
			emptydetector = np.empty((minydetectorsize,minxdetectorsize))




	 
		splitarrayx2 = np.array_split(emptydetector, cutoutsonyaxis, axis=0)


		#this section cuts off the last row to make all of the sections even. this means cutting off 1 bandwidth 
		#this goes down the rows 
		min_rows = min(arr.shape[0] for arr in splitarrayx2)
		min_cols = min(arr.shape[1] for arr in splitarrayx2)
		

		splitarrayx = []
		for section in splitarrayx2:
			
			#if the sections for each row  are not the same size, it cuts them to make sure they are 
			while section.shape[0] > min_rows:
				section = np.delete(section, section.shape[0] -1 , 0)
				

			while section.shape[1] > min_cols:
				section = np.delete(section, section.shape[1] -1 , 1)

			splitarrayx.append(section)
		
				


		splitdetector = []
		
		#this cuts every row into the correct number of columns
		for i in range(0, cutoutsonyaxis):

			splitarray = np.array_split(splitarrayx[i], cutoutsonxaxis, axis=1)



			
			
			splitdetector.append(
				splitarray
				
				)
		


		#earlier in the code we compared row by row to make sure they were totally the same size
		#this compares column to column across each section
		#doing both may be overkill, and we could potentially get away with just doing what's below, but why not do the overkill?

		step = 0
		
		
		for q in range(0, cutoutsonxaxis):
				

			for z in range(0, cutoutsonyaxis ):
				#goes through all cuts

				#print("section shape", splitdetector[z][q].shape, step)
				step +=1

				
				
				if q+1 >= cutoutsonxaxis:
					break
				#makes sure the sizes match for every section
				if splitdetector[z][q].shape != splitdetector[z][q+1].shape:
					print("Sizes for each detector do not match. Shape 1 = ", splitdetector[z][q].shape, 'Shape 2 =', splitdetector[z][q+1].shape, "Trimming shape 1 to be the same size..." )

				
					#the q+1 shape is likely to be smaller since it is on the end, that's how the array split works 
					shape_row = splitdetector[z][q+1].shape[0]
					shape_col = splitdetector[z][q+1].shape[1]

					#if row or column is not the smae, cuts a pixel off until it is
					#Note: "cutting off a pixel" in the detector is the same as just not using that space in real life
					#they may not be able to make a detector of that exact size, so just size up to the next biggest detector and don't use all the space
					#(or do, it's your detector after all im not the boss)
					while splitdetector[z][q].shape[0] > shape_row:

						splitdetector[z][q] = np.delete(splitdetector[z][q], splitdetector[z][q].shape[0] -1 , 0)
						


					while splitdetector[z][q].shape[1] > shape_col:
						splitdetector[z][q] = np.delete(splitdetector[z][q], splitdetector[z][q].shape[1] -1 , 1)

				#print("shape 1 update ", z, q, splitdetector[z][q].shape)

		

		

		
		
		yshape = splitdetector[0][0].shape
		#first sets your height, second sets your width
		


		#the -1 just makes sure that we aren't useing the last section
		#again, overkill since all of the sections should be the same but why not overkill?

		widthof1slice =  splitdetector[cutoutsonyaxis-1][cutoutsonxaxis-1].shape[1] #pixels

		
		heightof1slice = splitdetector[cutoutsonyaxis-1][cutoutsonxaxis-1].shape[0]
	

		#recalculating spatial and spectral resolution given cutouts
		#since the bandwidth actually goes over 2 slices, that's what we account for (nyquist sampling)
		detectorbandpass = bandwidth * (widthof1slice / 2)

		
		detectorFOV =  detectorpixelscale * heightof1slice


		print(' ')
		print(' ')
		print("Given your inputs: ")
		print("pixel pitch (microns): ", pixelpitch*1000)
		print("pixel scale (mas/pixel): ", detectorpixelscale/ 1e-3 )
		print("Desired FOV (arcsec squared): ", minFOV)
		print("Desired bandpass (nm): ", bandpass)
		print("R: ",R)
		print("Detector Size (mm): ", detectorsize.shape)
		print("Bandpass Start (nm): ", bandpassstart)
		print(' ')
		print("After setting up detector, these are the actual specs: ")
		print("Minimum detector size (mm) (this is what is used if you chose not to optimize): ", minydetectorsize*pixelpitch, ", " , minxdetectorsize*pixelpitch)
		print("Bandpass: ", detectorbandpass)
		print("FOV (arcsec squared): ", detectorFOV)
		print("total number of boxes (should be >= IFU slices): ", len(splitdetector)*len(splitarray))
		print("pixel pitch, x-axis (mas/pixel): ", (detectorFOV/numberIFUslices) * 1000)
		print("pixel scale, y-axis (mas/pixel): ",		(detectorFOV/heightof1slice) * 1000)
		print(' ')
		print(' ')



	time1 = datetime.now()	
	print("just got through input transform completely, current time:", time1.strftime("%H:%M:%S"))	
	return  detectorsizeinpixels, minxdetectorsize, minydetectorsize, bandwidth, pixelpitch, detectorpixelscale, R ,bandpass, possibledetector, lamb, bandpassstart, numberIFUslices, detectorbandpass, detectorFOV, splitdetector, cutoutsonxaxis, cutoutsonyaxis, widthof1slice

# [bytscl] -----------------------------------------------
#I actually don't think i use this function in the final version of PyISH
#but have fun with it if you want it
# ---------------------------------------------------------------------------

def bytscl_python(array, min_val=None, max_val=None, top=255):
    """
    Replicates the functionality of IDL's BYTSCL function in Python.

    Args:
        array (numpy.ndarray): The input array to be scaled.
        min_val (float, optional): The minimum value of the array to consider for scaling.
                                   If None, the actual minimum of the array is used.
        max_val (float, optional): The maximum value of the array to consider for scaling.
                                   If None, the actual maximum of the array is used.
        top (int, optional): The maximum value of the output byte range. Defaults to 255.

    Returns:
        numpy.ndarray: The scaled array as uint8 (byte) type.
    """
    if min_val is None:
        min_val = np.nanmin(array)  # Use nanmin to handle NaN values
    if max_val is None:
        max_val = np.nanmax(array)  # Use nanmax to handle NaN values

    # Handle cases where max_val equals min_val to avoid division by zero
    if max_val == min_val:
        return np.full_like(array, top // 2, dtype=np.uint8) # Return middle value if range is zero

    scaled_array = (array - min_val) * (top / (max_val - min_val))
    scaled_array = np.clip(scaled_array, 0, top) # Clip to ensure values are within 0 and top
    return scaled_array.astype(np.uint8)

def alogscale(image, minval, maxval, auto = 0):
    if auto == 1:
        med = np.nanmedian(image)
        sig = np.nanstd(image)
        maxval = (med + (10*sig)) < np.nanmax(image)
        minval = (med - (2*sig)) > np.nanmin(image)
    offset = minval - (maxval - minval) * .01
    scaled_image = bytscl_python(np.log10(image - offset), min_val = np.log10(minval - offset), max_val = np.log10(maxval-offset))
    print("min = " + str(np.log10(minval - offset)))
    print("max = " + str(np.log10(maxval-offset)))
    return(scaled_image)



# [Final Products] -----------------------------------------------
#THIS BAD BOY DOES A LOT OF WORK!
#this is the functuion that takes outputs from all previous functions and puts it together
#it takes the spectra from the input 3d cubes and places it in the detector sections we organized
#it then makes all the data products that are downloaded
# ---------------------------------------------------------------------------

def finalproducts(cutouts_dict, emission_wavelength_dict, lines, rebinned_cutouts_dict, detectorsizeinpixels, bandwidth,pixelpitch, detectorpixelscale, R,bandpass, possibledetector, lamb, bandpassstart, contcutouts, albedowavelength, albedoflux,detectorbandpass, detectorFOV, splitdetector, cutoutsonxaxis, cutoutsonyaxis, fulldetectorfits, threedcube, detectorsectionsfits):


	widthof1slice = splitdetector[0][0].shape[1] #pixels
	heightof1slice = splitdetector[0][0].shape[0]

	rebinnedcontinuumcutouts = [[] for _ in range(0,len(contcutouts))]

	sim_inputs = pd.read_excel('inputdataIFUSim.xlsx')
	#This is the height in arcsec of the hubble data. Look at word doc from sept 4 to see how I did this
	arcsecdata = sim_inputs.loc[0, 'arcsecspectra']

	inputs = pd.read_excel('InputInformationIFU.xlsx')

	

	darkrate = inputs.loc[0, 'dark rate']  #counts/sec/cm
	#pixelpitch in mm
	ExpTime = inputs.loc[0, 'ExpTime']

	pixelpitchcm = pixelpitch/10 #in cm
	darkratetomult = darkrate * (pixelpitchcm**2)*ExpTime #this is in counts
	


	pixelsasdetector = int(arcsecdata / detectorpixelscale)

	detector_cont = []	
	verticalstripsofdetector = []

	detectorsections = []

	

	

	

	#split detector is the array that has every detector section
	#i make a copy of it for the continuum data only
	#then add the emission lines to a seperate detector
	#then add the two together for one large detector
	continuum_detector_list = copy.deepcopy(splitdetector)




	#for the spectra from the continuum, we want to add spectra across the x-axis of the detector.
	#q and z go through every section in the detector (which corresponds to every IFS slice)
	slices = 0
	for q in range(0, cutoutsonxaxis):
				
		#if we reach a section that is not associated with a slice, break
		if slices > numberIFUslices:
			break

		for z in range(0, cutoutsonyaxis ):
			pixelplacement = 0
			
			#goes through every image in the input 3D cube, puts the desired slice for every image across bandpass
			for image in range(0, len(contcutouts)):
				
				
			
				#once slice within an image is reshaped, as each IFS slice is collapsed 
				cutout_continuum = np.sum(contcutouts[image][slices], axis = 1)

				cutout_continuum = np.reshape(cutout_continuum, (cutout_continuum.shape[0], 1))
				
			
				#duplicates to get every emission line 2 pixels wide. this is how nyquist sampling works 
				cutout_continuum = np.concatenate((cutout_continuum, cutout_continuum), axis = 1)
				
	
				
				#print(pixelplacement, pixelplacement+1)
				try:
					
					#we are placing the spectra from that slice in the correct pixel on the detector.
					np.put_along_axis(continuum_detector_list[z][q], np.array([[pixelplacement, pixelplacement + 1]]), [cutout_continuum], axis=1)
				except:
					print("continuum was not inserted at pixel ", pixelplacement,'/ ', (len(contcutouts) - 1)* 2, "likely because x-axis is an odd number of pixels and we are Nyquist Sampling." )
					print("Section (y, x): ", z, q)
					continue
			
				#goes every other since we are putting stuff across 2 pixels 
				#this is going through the x axis of every section on the detector
				pixelplacement = pixelplacement + 2

				rebinnedcontinuumcutouts[image].append(
					cutout_continuum)

			#pdb.set_trace()
		


					
			#need to match the IFS slice we are inputing to the detector section we are on	
			slices = slices + 1 
					
			if slices >= numberIFUslices:
						
				break

			
	
	#this is the loop through emission lines, to place them in the correct locations 

	for iteration, wavelengthlist, wavelengthimage in zip(range(0, len(emission_wavelength_dict)), emission_wavelength_dict , cutouts_dict.values()):

	



	  


		counter = 0

		

		p = 0
		for j in range(0, cutoutsonxaxis):
				

			if p > numberIFUslices-1:
					break

			for k in range(0, cutoutsonyaxis):
				counter = 0

				if p > numberIFUslices-1:
						
						break
				
				#Rather than going through the entire x axis of the detector section, as the continuum does, we systematically place the emission lines were they line up in the detector section
				for wavelength in wavelengthlist:

					

					
					#this is converting the wavelength for the emission line to the x-axis location of the detector section pixel that wavelength is associated with
					emissionlineplacement = int((wavelength- bandpassstart)/bandwidth) * 2 

					#used commented out stuff for debugging
					#print("ITERATION", iteration)
					#print("emissionlineplacement", emissionlineplacement)
					
					
				

					#same as continuum, collapse spectra to reflect how it is seen in the detecot
					rebinnedHRIcutout = np.sum(wavelengthimage[counter][p], axis = 1)
					rebinnedHRIcutout = np.reshape(rebinnedHRIcutout, (rebinnedHRIcutout.shape[0], 1))

					

			

					#duplicate such that each emission line goes over 2 pixels
					rebinnedHRIcutout = np.concatenate((rebinnedHRIcutout, rebinnedHRIcutout), axis = 1)

					#add 0's until it matches height dimensio of detector
					#rebinnedHRIcutout  = rebin(rebinnedHRIcutout, splitdetector[k][j].shape[0])
					#place in detector
					np.put_along_axis(splitdetector[k][j], np.array([[emissionlineplacement, emissionlineplacement + 1]]), [rebinnedHRIcutout], axis=1)

					counter = counter + 1
		
					
				#This just makes sure I am not going over the amount of slices I have 
				p = p + 1 
					
				if p > numberIFUslices-1:
						
					break
				if iteration == len(emission_wavelength_dict) -1:

					#when we get to the last emission line in the dictionary, we add our emission + continuum detectors together to get one detector
					#we also add in dark rate and noise
					
					splitdetector[k][j] = np.add(splitdetector[k][j], continuum_detector_list[k][j])
						

					splitdetector[k][j] = splitdetector[k][j] + darkratetomult
					
					splitdetector[k][j] = np.random.poisson(splitdetector[k][j])

					#np.floor insures that all values are integers (since they are photons)
					splitdetector[k][j] = np.floor(splitdetector[k][j])
					
					
					#save them
					detectorsections.append(splitdetector[k][j])







	if detectorsectionsfits == True:

		#this will save a fits file that has the axes (detector x axis, detector y axis, section)

		hdudetector = fits.PrimaryHDU(detectorsections)
		hduld = fits.HDUList([hdudetector])
		hduld.writeto("detectorsections.fits", overwrite=True)
	















	#this goes through every detector section (aka split detector), which now has data in it, and stitches them together to form one large full format detector
	#the 50 is the width (in pixels) of the dividers seperating each section, and 15 is the value of those dividers
	# if your detector has spectra that contain photons MUCH greater than 15, increase 15 to some number so you can see the section dividers
	#or, uncomment the things I have commented out within the j and k loop to just not include dividers
	
	#this first loop stitches all of the vertical sections of the detector together (all of the columns)	
	for j in range(0, cutoutsonxaxis):

		for k in range(0, cutoutsonyaxis ):

			
			
							
			if k == 0:

				
							
				addedtogether =  np.concatenate((splitdetector[k][j], np.full((50,splitdetector[k][j].shape[1]), 15)), axis = 0)
				addedtogether =  np.concatenate((addedtogether, splitdetector[k + 1][j]), axis = 0)
				#uncomment this below, and comment out 2 lines above to not include the dividers between sections
				#addedtogether =  np.concatenate((splitdetector[k][j], splitdetector[k + 1][j]), axis = 0)
							
			else: 
							
				while k + 1 < cutoutsonyaxis:

					addedtogether =  np.concatenate((addedtogether, np.full((50,splitdetector[k][j].shape[1]), 15)), axis = 0)
					addedtogether =  np.concatenate((addedtogether, splitdetector[k + 1][j]), axis = 0)
					#uncomment this below, and comment out 2 lines above to not include the dividers between sections
					#addedtogether =  np.concatenate((addedtogether, splitdetector[k + 1][j]), axis = 0)
						
								
								

					break

				
		#saves the vertical strips to stitch back together fully after
		verticalstripsofdetector.append(
			addedtogether)


		#This group ensures that the heights for each detector are the same
		#if they are not, means that the amount of cutouts is uneven, so add an empty section of the detector
	for l in range (1,len(verticalstripsofdetector)):

		splitdetector[cutoutsonyaxis-1][cutoutsonxaxis-1].fill(0)
			
		while verticalstripsofdetector[l].shape[0] != verticalstripsofdetector[l-1].shape[0]:
			verticalstripsofdetector[l] = np.concatenate((verticalstripsofdetector[l], np.full((50,splitdetector[k][j].shape[1]), 15)), axis = 0)
			verticalstripsofdetector[l] = np.concatenate((verticalstripsofdetector[l], splitdetector[cutoutsonyaxis-1][cutoutsonxaxis-1]), axis= 0)
			



	#This section adds those columns together into the full detector
	#again, uncomment the code in the loop, and comment out the stuff already there to stitch the detector together without those section dividers
	for l in range (0,len(verticalstripsofdetector)):

		if l == 0:

			finaldetector =  np.concatenate((verticalstripsofdetector[l], np.full((verticalstripsofdetector[l].shape[0],50), 15)), axis = 1)
			finaldetector =  np.concatenate((finaldetector, verticalstripsofdetector[l+1]), axis = 1)
			#finaldetector =  np.concatenate((verticalstripsofdetector[l], verticalstripsofdetector[l+1]), axis = 1)
			
		else: 
						
			while l + 1 < cutoutsonxaxis:
							
				finaldetector =  np.concatenate((finaldetector, np.full((verticalstripsofdetector[l].shape[0],50), 15)), axis = 1)
				finaldetector =  np.concatenate((finaldetector, verticalstripsofdetector[l+1]), axis = 1)

				#finaldetector =  np.concatenate((finaldetector, verticalstripsofdetector[l+1]), axis = 1)
				
				break

	

	#print("final detector shape", finaldetector.shape)

	#below prints these out to the user
	'''
	print("here we can print the detector sections")
	finaldetector[finaldetector != 15] = 0
	plt.imshow(finaldetector)
	plt.xticks([])
	plt.yticks([])
	plt.xticks([])

	#plt.savefig('/Users/gsweetak/Documents/JATIS2025/revisions/detectorsectionsafter.png', dpi = 600)
	
	plt.show()

	pdb.set_trace()
	'''

	#this is what saves the full detector as a fits file
	if fulldetectorfits == True:
	
		hdu = fits.PrimaryHDU(finaldetector)
		hdul = fits.HDUList([hdu])
		hdul.writeto("stitcheddetector.fits", overwrite=True) 




#Stitching 3D cube together

	endlength = int(detectorsections[0].shape[1]/2)

	stripsofslice = [[] for _ in range(0,endlength)]
	slices1 = 0
	for section in detectorsections:
			pixel_section = 0
			step1 = 0
			if slices1 > numberIFUslices:
				break
			
			
			for step in range(0, section.shape[1], 2):

			#go through x axis of each detector section, and add them into a list
			#this repeats for every section, which essentially reccreates the input images as output images (from what the detecor sees)
			#so each [step] in this stripsofslice list corresponds to a wavelength, and the information in each [step] forms the spatial maps
				strip = section
				strip = strip[:, step:step + 1]
				
				

				stripsofslice[step1].append(
					strip)

				step1 = step1 + 1

				if step1 >= endlength:
					break
					
				
			slices1 = slices1 + 1 

			
					
			if slices1 >= numberIFUslices:		
				break
	restitchedimages = []

	#this ensures everything is the same size (overkill)
	for v in range(0, len(stripsofslice)):
		for individual_slice in range(1, len(stripsofslice[v])):



			if stripsofslice[v][individual_slice].shape[0] != stripsofslice[v][individual_slice-1].shape[0]:
					stripsofslice[v][individual_slice] = rebin(stripsofslice[v][individual_slice], stripsofslice[v][individual_slice-1].shape[0] )
					
		#concatanating them is just forming each stripsofslice[v] into one spatial map per v, rather than a bunch of arrays per v
		twodmap = np.concatenate(stripsofslice[v], axis = 1)

		#going to attempt resizing so that the cube is a box
		#first reshapes to pixels as detector wide
		#this is changing the units to (x,y) = (pixels, pixels) rather than (numberofifsslices, pixels)
		twodmap = resample(twodmap, (twodmap.shape[0], heightof1slice), method = 'linear')
		twodmap = np.clip(twodmap, 0, None)
		twodmap = np.floor(twodmap)


			

		restitchedimages.append(twodmap)

	#saves the 3d cube
	if threedcube == True:

		hdu2 = fits.PrimaryHDU(restitchedimages)
		hdul2 = fits.HDUList([hdu2])
		hdul2.writeto("3dcube_HWOsim.fits", overwrite=True) 




	time1 = datetime.now()	
	print("just got through final products completely, current time:", time1.strftime("%H:%M:%S"))	

	#this goes through picking out a particular spaxel and showing a spectrum of that spaxel 

	showing1dspectrum = input("Would you like to see a one dimensional spectrum of a particular spaxal in the cube? If so, type y. Otherwise, type n.")
	if showing1dspectrum == "y":
		#shows the user a spatial map of one emission line and prompt to pick a coordinate from that map
		plt.imshow(restitchedimages[int(emissionlineplacement/2)])
		plt.colorbar()
		plt.xticks([])
		plt.yticks([])
		plt.show(block=False)
		plt.savefig('europamap.png', dpi = 600)

		print("What spaxal's 1D spectrum would you like to see? Enter the x and y coordinates below:")
		x = int(input("x: "))
		y = int(input("y: "))

		plt.close()
		onedspectrum = []

		#rhis is going through each spatial map in the 3D cube 'restitched images' and extracting the flux at the requested pixel
		for spectralmap in range(0, len(restitchedimages)):

			fluxval = restitchedimages[spectralmap]
			fluxval = float(fluxval[y][x])
			onedspectrum.append(fluxval)
		#pdb.set_trace()

		#plots the flux
		plt.plot(np.arange(bandpassstart + bandwidth, bandpassstart + detectorbandpass, bandwidth), onedspectrum)
		plt.xlabel('Wavelength [nm]')
		plt.ylabel("Photons/nm")
		#plt.grid()
		#plt.yscale('log')
		plt.savefig('1dspectrum.png', dpi = 600)
		plt.show()
		pdb.set_trace()
		plt.close()



	






# ==============================================================================
#                        M A I N   F U N C T I O N  
# ==============================================================================
# Everything below this block contains the core execution logic of the script.
# I know the variables are disgusting looking, I apologize
# she's ugly but she works
# ==============================================================================





#start by using function to understand if user inputs are possible to accomadate with an IFS
detectorsizeinpixels, minxdetectorsize, minydetectorsize, bandwidth, pixelpitch, detectorpixelscale, R ,bandpass, possibledetector, lamb, bandpassstart, numberIFUslices, detectorbandpass, detectorFOV, splitdetector, cutoutsonxaxis, cutoutsonyaxis, widthof1slice= inputtransform()


#if they are not possible, can't do anything else
if possibledetector == 0:
	print("This detector setup does not work, so I cannot continue.")
else:

	#if they are possible, go through functions earlier
	data_list, albedowavelength, albedoflux, contdata_listR, emission_wavelength_dict, widthemissionlines, emissionline_wl, lines= imageprocessing(widthof1slice, bandwidth, detectorbandpass, bandpassstart, detectorFOV)

	scaledHRIdata, contdata_listphoton= telescopecalibration(data_list, emission_wavelength_dict, contdata_listR,  bandpassstart, detectorbandpass, bandwidth, widthemissionlines, emissionline_wl, detectorpixelscale, splitdetector, cutoutsonxaxis, cutoutsonyaxis)

	contcutouts = continuumslicer(contdata_listphoton, numberIFUslices)

	cutouts_dict, rebinned_cutouts_dict= ifuslicer(scaledHRIdata, numberIFUslices, lines)

	
	#--------------IMPORTANT NOTE------------
	#fulldetectorfits = True, means save a .fits file of your final IFS detector filled with spectra (3D)
	#threedcube = True, hopefully the variable name is self explanatory. (x,y,z) = (pixels, pixels, wavelength (pixel units though, must do wavelength solution yourself))
	#detectorsectionfits = True, a stitched .fits file of every section in your detector (2D)
	#switch these to False if you don't want to save them!

	finalproducts(cutouts_dict,

	emission_wavelength_dict,
	lines,
	rebinned_cutouts_dict,
	 detectorsizeinpixels,
	  bandwidth,
	  pixelpitch, 
	  detectorpixelscale,
	   R,
	   bandpass,
	    possibledetector, 
	    lamb,
	     bandpassstart, 
	     contcutouts,
	      albedowavelength, 
	      albedoflux, 
	      detectorbandpass, 
	      detectorFOV, 
	      splitdetector, 
	      cutoutsonxaxis, 
	      cutoutsonyaxis, 
	      
	      fulldetectorfits = True, 
	      threedcube = True, 
	      detectorsectionsfits = True
	  )