# -*- coding: utf-8 -*-
"""
Created on Mon Apr 27 12:11:17 2026

@author: Alberto
"""


# =============================================================================
# 1 - Import libraries and set-up graphics style
# =============================================================================


import matplotlib.pyplot as plt #Used Matplotlib version 3.10.8
import numpy as np #Used numpy version 1.26.4

import pandas as pd #Used pandas version 3.0.2

import seaborn as sns #Used seaborn version 0.13.2
import os

from scipy.optimize import curve_fit #Used scipy version 1.11.1
from scipy.special import erf

sns.set_theme()
sns.set_style("ticks")
sns.set_context("paper",font_scale=1.5)
sns.set_palette("tab10")

# =============================================================================
# 2 - Define working directory, fitting function (Solution to 1D diffusion equation) and create folders for saving data
# =============================================================================

os.chdir(r"G:\OneDrive Doctorado\ALBA\ID2024098669\Analysis\Maps one step\AuNPs_22nm_PolyT22_PBS_11.52nM")

def Diffusion_erf_2(x, D, t0,C_0,x_0):
    """
    Solution to the 1D Diffusion Equation with constant initial concentration over time for fitting

    Parameters
    ----------
    x : Distance to initial point
    D : Diffusion Coefficient
    t0 : Initial Diffusion time (to account for experimental time delays running the experiment)
    C_0 : Initial concentration (Concentration at the NP channel)
    x_0 : NP-Matrigel frontier position (to account for spatial displacement)

    Returns
    -------
    y : Concentration profile

    """
    
    y = C_0/2*(1+erf((-x+x_0)/np.sqrt(4*D*(t-t0))))
    
    
    return y

[os.makedirs(folder,exist_ok=True) for folder in ["Results","Image Sequence/Au","Image Sequence/Fe"]]



# =============================================================================
# 3 - Initialize variables 
# =============================================================================

t = 0 #Time

img_Au = np.zeros((7,51,51)) #Time-lapse image
img_Fe = np.zeros((7,51,51)) #Time-lapse image


# =============================================================================
# 4 - Read data files, create and save fluorescence maps
# =============================================================================


files = os.listdir("Experimental Data")

files.sort()

print(files)
for file in files:
    
    #Read data file
    map_data = pd.read_table(f"Experimental Data/{file}", skiprows = 33,skipfooter=1,sep="\s+",engine="python") 
    
    #Rename columns
    map_data = map_data.set_axis(['Pt_No', 'dif_ms', 'mtp_sx', 'energyc', 'Iring', 'x_ch1_roi1', 'x_ch2_roi1', 'x_ch3_roi1', 'x_ch4_roi1', 'x_ch5_roi1', 'x_ch6_roi1',  'x_ch7_roi1', 'x_ch8_roi1', 'x_ch9_roi1', 'x_ch10_roi1', 'x_ch11_roi1',       'x_ch12_roi1', 'x_ch13_roi1', 'x_ch1_roi2', 'x_ch2_roi2', 'x_ch3_roi2',       'x_ch4_roi2', 'x_ch5_roi2', 'x_ch6_roi2', 'x_ch7_roi2', 'x_ch8_roi2',       'x_ch9_roi2', 'x_ch10_roi2', 'x_ch11_roi2', 'x_ch12_roi2',       'x_ch13_roi2', 'mtp_ic_timer', 'mtp_ic0', 'mtp_ic1', 'dt',"1"],axis=1)
    
    #Get motor positions
    x_pos = set(map_data.loc[:,"mtp_sx"])
    y_pos = set(map_data.loc[:,"dif_ms"])
    
    #Sum fluorescence intensity from all 13 SDD channels
    map_data.loc[:,"sum_roi_1"] = map_data.loc[:,"x_ch1_roi1"]
    map_data.loc[:,"sum_roi_2"] = map_data.loc[:,"x_ch1_roi2"]
    
    for i in range(2,14):
        
        map_data.loc[:,"sum_roi_1"] += map_data.loc[:,f"x_ch{i}_roi1"]
        map_data.loc[:,"sum_roi_2"] += map_data.loc[:,f"x_ch{i}_roi2"]
    
    #Create map
    data_roi1 = np.zeros((len(x_pos),len(y_pos)))
    data_roi2 = np.zeros((len(x_pos),len(y_pos)))
    
    
    x_pos_sorted = list(x_pos)
    x_pos_sorted.sort()
    
    
    for i,x in enumerate(x_pos_sorted):
        
        
        
        data_roi1[:,i] = map_data.loc[map_data.loc[:,"mtp_sx"]==x,"sum_roi_1"]
    
        data_roi2[:,i] = map_data.loc[map_data.loc[:,"mtp_sx"]==x,"sum_roi_2"]
        
    

    #Plot map
    fig1, (ax1,ax2) = plt.subplots(1,2)
    
    
    ax1.imshow(data_roi1,cmap="jet", vmin = 0.9*np.min(data_roi1), vmax = 1.2*np.max(data_roi1))
    
    ax1.set_title("Fe")
    
    ax2.imshow(data_roi2,cmap="jet", vmin = 0.9*np.min(data_roi2),vmax = 1.2*np.max(data_roi2))
    ax2.set_title("Au")
    
    plt.suptitle(f"t = {t} mins")
    
    plt.show()
    plt.close()
    
    #Save map data to txt file (Can be opened in fiji/ImageJ to convert it to tiff)
    save_map_data_Fe = pd.DataFrame(data_roi1)
    save_map_data_Au = pd.DataFrame(data_roi2)
    
    
    save_map_data_Fe.to_csv(f"Image Sequence/Fe/Fe t = {t} mins.txt",sep="\t",header=None,index=None)
    save_map_data_Au.to_csv(f"Image Sequence/Au/Au t = {t} mins.txt",sep="\t",header=None,index=None)
    

    #Append map data to time-lapse image
    img_Au[int(t/20),:,:] = data_roi2
    img_Fe[int(t/20),:,:] = data_roi1

    t+=20
    
    print(file)



# =============================================================================
# 5 - Begin image analysis
# =============================================================================

#Choose signal to use for analysis
element = input("Which element do you want to analyze? (Fe / Au): ")


if element == "Fe":
    img = img_Fe
elif element == "Au":
    img = img_Au

else:
    raise TypeError("Wrong element selected, only Fe or Au allowed. Maybe you misspelled, beware of Capital letters.")

#Create folder for saving analyis data
[os.makedirs(folder,exist_ok=True) for folder in [f"Plots/{element}/Maps",f"Plots/{element}/Parameters"]]


#Define line to obtain profile concentrations
cont=0

while cont != "y" :
    
    
    x = int(input("First Point X coord: "))
    
    y = int(input("First Point Y coord: "))
    
    yarray = np.arange(0,y+1,1)
    
    xarray = np.zeros(len(yarray))
    
    xarray+=x
    
    for i in range(img.shape[0]):
            
        plt.imshow(img[i,:,:],cmap="jet",vmin = 0.9*np.min(img[i,:,:]), vmax = 1.2*np.max(img[i,:,:]))
        plt.plot(xarray,yarray,"r")
        
        plt.show()
        plt.close()
    
    cont =input("Continue? (y/n): ")



#Initialize list to save fitted parameters 
Ds = []

C_0s = []

t_0s = []

x_0s = []

t = 0


colors = sns.color_palette("bright")
colors.as_hex()

fig1 = plt.figure(1)
ax1 = fig1.gca()




for i in range(img.shape[0]):
        
    fig2 = plt.figure(2)
    ax2 = fig2.gca()
    
    #Obtain Concentration profiles
    X = np.arange(0,y+1,1) #X-positions in pixels 1 pixel  == 44 um
    C = img[i,0:y+1,x]
    
    
    C_list = list(C)
    C_list.reverse()
    C = np.array(C_list)
    
    #Subtract background
    C-=np.min(C)
    
    #Normalize concentration
    C_norm = C/np.max(C)
    
    #Plot Concnetration Profiles
    if i !=0:
        ax1.plot(X*44,C_norm,"o",label = f"t = {t/60} min",color=colors[i],markerfacecolor="None")
        ax2.plot(X*44,C_norm,"o",label = f"t = {t/60} min",color=colors[i],markerfacecolor="None")
    
    #Perform fitting
    pars,cov = curve_fit(Diffusion_erf_2,X*44,C_norm,p0=[1,-30,1,1],bounds = ((0,-np.inf,0,-np.inf),(np.inf,0,np.inf,np.inf)))
    
    #Plot fitting
    X_fit = np.linspace(X[0],X[-1],10000)
    if i !=0:
        ax1.plot(X_fit*44,Diffusion_erf_2(X_fit*44,*pars),"--",color=colors[i])
        ax2.plot(X_fit*44,Diffusion_erf_2(X_fit*44,*pars),"--",color=colors[i])
    
    #Extract fitted parameters and append to saving lists
    D,t0,C_0,x_0 = pars
    
   
    Ds.append(D)
    
    C_0s.append(C_0)
    
    t_0s.append(t0)
    
    x_0s.append(x_0)
           

        
    ax2.set_xlabel("Distance ($\mu$m)")
    ax2.set_ylabel("Normalized Intensity")
    
    ax2.set_title(f"t = {int(t/60)} min")
    
    #Save concentration profile and fit plot for each time point
    
    fig2.savefig(f"Plots/{element}/t = {int(t/60)} min.tif",dpi=300,bbox_inches="tight")
    fig2.show()
    plt.close(fig2)

    t+=20*60

    

ax1.set_xlabel("Distance ($\mu$m)")
ax1.set_ylabel("Normalized Intensity")
ax1.legend()
ax1.set_xlim([-61,1300])

#Save all concentration profiles and fits in one plot
fig1.savefig(f"Plots/{element}/All.tif",dpi=300,bbox_inches="tight")
fig1.show()
plt.close(fig1)


    
print(Ds,t_0s,C_0s,x_0s)


#Plot fitting parameters vs time
times = np.linspace(0,120,7)

plt.plot(times,C_0s,".-")
plt.title("Fitting Parameter C$_0$")
plt.ylabel("C_$0$")
plt.xlabel("Time (min)")
plt.savefig(f"Plots/{element}/Parameters/C_0s.tif",bbox_inches="tight")
plt.show()
plt.close()


plt.plot(times,t_0s,".-")
plt.title("Fitting Parameter t$_0$")
plt.ylabel("t_$0$ (s)")
plt.xlabel("Time (min)")
plt.savefig(f"Plots/{element}/Parameters/t_0s.tif",bbox_inches="tight")
plt.show()
plt.close()


plt.plot(times,x_0s,".-")
plt.title("Fitting Parameter x$_0$")
plt.ylabel("x_$0$ ($\mu$m)")
plt.xlabel("Time (min)")
plt.savefig(f"Plots/{element}/Parameters/x_0s.tif",bbox_inches="tight")
plt.show()
plt.close()



plt.plot(times,Ds,".-")
plt.title("Fitting Parameter D")
plt.xlabel("Time (min)")
plt.ylabel("D ($\mu$m$^2$/s)")
plt.savefig(f"Plots/{element}/Parameters/Ds.tif",bbox_inches="tight")
plt.show()
plt.close()


#Save fitting results as txt file

log = pd.DataFrame()

log.loc[0,"Parameter"] = "First Point Coordinates (pixels x pixels):"

log.loc[0,"Value"] = f"({x},{y})"

log.loc[1,"Parameter"] = "Average D ($mu$m$^2$/s)"

log.loc[1,"Value"] = f"{np.mean(Ds[1:])} $\pm$ {np.std(Ds[1:])}"


log.loc[2,"Parameter"] = "Average C_0"

log.loc[2,"Value"] = f"{np.mean(C_0s[1:])} $\pm$ {np.std(C_0s[1:])}"


log.loc[3,"Parameter"] = "Average t_0 (s)"

log.loc[3,"Value"] = f"{np.mean(t_0s[1:])} $\pm$ {np.std(t_0s[1:])}"


log.loc[4,"Parameter"] = "Average x_0 ($mu$m)"

log.loc[4,"Value"] = f"{np.mean(x_0s[1:])} $\pm$ {np.std(x_0s[1:])}"


log["Value"] = log["Value"].astype('object')

log.at[5,"Parameter"] ="Ds"

log.at[5,"Value"] = Ds[1:]


log.at[6,"Parameter"] ="x_0s"

log.at[6,"Value"] = x_0s[1:]


log.at[7,"Parameter"] ="t_0s"

log.at[7,"Value"] = t_0s[1:]


log.at[8,"Parameter"] ="C_0s"

log.at[8,"Value"] = C_0s[1:]


print(log)

log.to_csv(f"Results/Results_{element}.txt",sep="\t",index=None)




