#!/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),
                         ('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,mua,mus,g):
    StructMilieu['Epaisseur']=epaisseur
    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,dz=fonctionPhase(StructMilieu,Structphoton)
    Structphoton['Direction'][0,0]=dx
    Structphoton['Direction'][0,1]=dz
    Structphoton['Poids']=Structphoton['Poids']-(mua/(mua+mus))*Structphoton['Poids']
    Structphoton=mortPhoton(Structphoton,StructMilieu)
    return Structphoton

def viePhoton(Structphoton,StructMilieu):
    TrajX=[]  # init liste des positions successives x du photon
    TrajZ=[]  # init liste des positions successives z du photon
    PoidsListe=[] # init liste des poids w du photon
    TrajX.append(float(Structphoton['Position'][0,0])) # ajout du x=0 initial
    TrajZ.append(float(Structphoton['Position'][0,1])) # ajout du z=0 initial
    PoidsListe.append(Structphoton['Poids'][0]) # ajout du w=1 initial
    ## compléter ici
    
    TrajListe=[TrajX,TrajZ]
    return TrajListe,PoidsListe

def reflectance_transmittance(StructSimu,StructMilieu):
    # compléter ici
    return StructSimu




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