#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Oct  6 15:21:53 2023

@author: victorcolas
"""

import numpy as np
import random as rd
import matplotlib.pyplot as plt


Structphoton = np.zeros(1, dtype = [('Alive', bool),
                         ('Position', float, (2, 1)),
                         ('Direction', float, (2, 1)),('Poids',float)])


StructMilieu = np.zeros(1, dtype = [('Epaisseur', float),
                         ('Indice', float),
                         ('Mua', float),('Mus',float),('G',float)])

StructSimu = np.zeros(1, dtype = [('NbPhoton', int),
                         ('R', float),
                         ('T', float)])


def initStrucPhoton(Structphoton):
    Structphoton['Alive']=True
    Structphoton['Position']=[[0],[0]]
    Structphoton['Direction']=[[0],[1]]
    Structphoton['Poids']=1
    return Structphoton
    
def initMilieu(StructMilieu,epaisseur,indice,mua,mus,g):
    StructMilieu['Epaisseur']=epaisseur
    StructMilieu['Indice']=indice
    StructMilieu['Mua']=mua
    StructMilieu['Mus']=mus
    StructMilieu['G']=g
    return StructMilieu

def initSimu(StructSimu,NbPhotons):
    StructSimu['NbPhoton']=NbPhotons
    StructSimu['R']=0
    StructSimu['T']=0
    return StructSimu

def fonctionPhase(StructMilieu,Structphoton):
    dxAncien=float(Structphoton['Direction'][0,0])
    dzAncien=float(Structphoton['Direction'][0,1])
    thetaAncien=np.arctan(dxAncien/dzAncien)
    g=float(StructMilieu['G'])
    u=rd.random()
    if g==0:
        Deviation=2*np.pi*u
        CosDeviation=np.cos(Deviation)
    else :
        temp=(1-g**2)/(1-g+2*g*u)
        CosDeviation=(1+g*g-temp*temp)/(2*g)
        Deviation=np.arccos(CosDeviation)
    theta=thetaAncien+Deviation
    dz=np.cos(theta)
    usigne=rd.random()
    if usigne<=0.5:
        signe=-1
    else:
        signe=1
    dx=signe*np.sqrt(1-dz**2)
    return dx,dz

def mortPhoton(Structphoton,StructMilieu):
    poids=float(Structphoton['Poids'])
    pos=Structphoton['Position']
    pos=pos[0]
    epaisseur=StructMilieu['Epaisseur'][0]
    z=(pos[1][0])
    if poids<0.01:
        Structphoton['Alive']=False
    elif z>epaisseur:
        Structphoton['Alive']=False
    elif z<0.:
        Structphoton['Alive']=False
    return Structphoton

def diffusion(Structphoton,StructMilieu):
    mus=StructMilieu['Mus']
    mua=StructMilieu['Mua']
    l=rd.expovariate(mus)
    Structphoton['Position']=Structphoton['Position']+l*Structphoton['Direction']
    dx,dy=fonctionPhase(StructMilieu,Structphoton)
    Structphoton['Direction'][0,0]=dx
    Structphoton['Direction'][0,1]=dy
    Structphoton['Poids']=Structphoton['Poids']-(mua/(mua+mus))*Structphoton['Poids']
    Structphoton=mortPhoton(Structphoton,StructMilieu)
    return Structphoton

def viePhoton(Structphoton,StructMilieu):
    TrajX=[]
    TrajZ=[]
    PoidsListe=[]
    TrajX.append(float(Structphoton['Position'][0,0]))
    TrajZ.append(float(Structphoton['Position'][0,1]))
    PoidsListe.append(Structphoton['Poids'][0])
    while Structphoton['Alive']==True :
        Structphoton=diffusion(Structphoton,StructMilieu)
        TrajX.append(float(Structphoton['Position'][0,0]))
        TrajZ.append(float(Structphoton['Position'][0,1]))
        PoidsListe.append(Structphoton['Poids'][0])
    
    TrajListe=[TrajX,TrajZ]
    return TrajListe,PoidsListe

def reflectance_transmittance(StructSimu,StructMilieu):
    nPhot=StructSimu['NbPhoton'][0]
    for k in range(nPhot):
        Structph=initStrucPhoton(Structphoton)
        TrajListe,PoidsListe=viePhoton(Structph,StructMilieu)
        if TrajListe[1][-1]<0:
            StructSimu['R']=StructSimu['R']+PoidsListe[-1]
        elif TrajListe[1][-1]>StructMilieu['Epaisseur'][0]: 
            StructSimu['T']=StructSimu['T']+PoidsListe[-1]
    StructSimu['R']=StructSimu['R']/nPhot
    StructSimu['T']=StructSimu['T']/nPhot
    return StructSimu




### Simu puis Résultats
# param simu
nPhot=10000
StructSimu=initSimu(StructSimu, nPhot)
# def des propriétés du milieu
epaisseur=1    #mm
n=1.4
mua=0.2          #/mm
mus=5         #/mm
g= 0.4
StructMilieu=initMilieu(StructMilieu, epaisseur,n,mua,mus,g)


### test des fonction
Structphoton=initStrucPhoton(Structphoton)
dx,dz=fonctionPhase(StructMilieu,Structphoton)


### Affichage traj photon
StructSimu=reflectance_transmittance(StructSimu,StructMilieu)
print(StructSimu['R'])
print(StructSimu['T'])
plt.figure()
for k in range(10):
    Structphoton=initStrucPhoton(Structphoton)
    TrajListe,PoidsListe=viePhoton(Structphoton,StructMilieu)
    plt.plot(TrajListe[1],TrajListe[0])
    
plt.grid(1)
plt.xlabel('z (mm)')
plt.ylabel('x (mm)')

# # R et T fonction de l'épaisseur
plt.figure()
nEp=20
epaiseurListe=np.linspace(0.1,3,nEp)
listeR=[]
listeT=[]
for k in range(nEp):
    print(k)
    StructMilieu=initMilieu(StructMilieu, epaiseurListe[k],n,mua,mus,g)
    StructSimu=reflectance_transmittance(StructSimu,StructMilieu)
    listeR.append(StructSimu['R'][0])
    listeT.append(StructSimu['T'][0])
  
plt.plot(epaiseurListe,listeR,label='R')
plt.plot(epaiseurListe,listeT,label='T')
plt.grid(1)
plt.legend()
plt.xlabel('Epaisseur (mm)')
plt.ylabel('R & T ')