transmission_modelling.py

# -*- coding: utf-8 -*-
"""Transmission modelling.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1Qorq5EQU91qs1znGXhUymkQPBuKDXqUm
"""

# install Pint if necessary

try:
    import pint
except ImportError:
    !pip install pint

# download modsim.py if necessary

from os.path import exists

filename = 'modsim.py'
if not exists(filename):
    from urllib.request import urlretrieve
    url = 'https://raw.githubusercontent.com/AllenDowney/ModSim/main/'
    local, _ = urlretrieve(url+filename, filename)
    print('Downloaded ' + local)

# import functions from modsim

from modsim import *
from numpy import *
import math
#Use to import pandas
import pandas as pd
#Use to import the file into google Colab drive
from google.colab import files 
#Use to import io, which opens the file from the Colab drive
import io
from google.colab import drive
uploaded = files.upload()
uploaded_file_name = list(uploaded.keys())[0]
df = pd.read_excel(io.BytesIO(uploaded.get(uploaded_file_name)))

uploaded = files.upload()
uploaded_file_name = list(uploaded.keys())[0]
uf = pd.read_excel(io.BytesIO(uploaded.get(uploaded_file_name)))

"""" Calculates the transmission efficiency using a gaussian model
   Arguments:
      theta: an angle in degrees representing the horizontal angle
      phi: angle in degrees representing the vertical angle
      sigma_L, sigma_R, sigma_B, sigma_T: parameters for the modified gaussian model
       chosen based on the dT
    Returns:  
      Transmission efficiency using the modified gaussian equation 
      for the given parameters """ 
def gaussian(theta, phi, sigma_L, sigma_R, sigma_B, sigma_T):
  ##from Jay models powerpoint
  if theta>=0 and phi>=0:
    return math.exp(-theta**2/(2*sigma_R**2))*math.exp(-phi**2/(2*sigma_T**2))
  elif theta>=0 and phi<=0:
    return math.exp(-theta**2/(2*sigma_R**2))*math.exp(-phi**2/(2*sigma_B**2))
  elif theta<=0 and phi<=0:
    return math.exp(-theta**2/(2*sigma_L**2))*math.exp(-phi**2/(2*sigma_B**2))
  else:
    return math.exp(-theta**2/(2*sigma_L**2))*math.exp(-phi**2/(2*sigma_T**2))

"""" Calculates the transmission efficiency using a modified gaussian model
   Arguments:
      theta: an angle in degrees representing the horizontal angle
      phi: angle in degrees representing the vertical angle
      sigma_L, sigma_R, sigma_B, sigma_T, theta0, N: parameters for the modified gaussian model
       chosen based on the dT
    Returns:  
      Transmission efficiency using the modified gaussian equation 
      for the given parameters """ 
def mod_gaussian(theta, phi, sigma_L, sigma_R, sigma_B, sigma_T, theta0, N):
  ##from Jay models powerpoint
  if theta>=theta0 and phi>=0:
    return N*math.exp(-(theta-theta0)**2/(2*sigma_R**2))*math.exp(-phi**2/(2*sigma_T**2))
  elif theta>=theta0 and phi<=0:
    return N*math.exp(-(theta-theta0)**2/(2*sigma_R**2))*math.exp(-phi**2/(2*sigma_B**2))
  elif theta<=theta0 and phi<=0:
    return N*math.exp(-(theta-theta0)**2/(2*sigma_L**2))*math.exp(-phi**2/(2*sigma_B**2))
  else:
    return N*math.exp(-(theta-theta0)**2/(2*sigma_L**2))*math.exp(-phi**2/(2*sigma_T**2))

"""" Calculates the transmission efficiency using a modified gaussian fixed sigma L model
   Arguments:
      theta: an angle in degrees representing the horizontal angle
      phi: angle in degrees representing the vertical angle
      sigma_R, sigma_B, sigma_T, theta0: parameters for the modified gaussian fixed sigma L model
       chosen based on the dT
    Returns:  
      Transmission efficiency using the modified gaussian fixed sigma L equation 
      for the given parameters """ 
def mod_gaussian_fixed_sigmaL(theta, phi, sigma_R, sigma_B, sigma_T, theta0):
  ##from Jay models powerpoint
  sigma_L = 10
  if theta>=theta0 and phi>=0:
    return math.exp(-(theta-theta0)**2/(2*sigma_R**2))*math.exp(-phi**2/(2*sigma_T**2))
  elif theta>=theta0 and phi<=0:
    return math.exp(-(theta-theta0)**2/(2*sigma_R**2))*math.exp(-phi**2/(2*sigma_B**2))
  elif theta<=theta0 and phi<=0:
    return math.exp(-(theta-theta0)**2/(2*sigma_L**2))*math.exp(-phi**2/(2*sigma_B**2))
  else:
    return math.exp(-(theta-theta0)**2/(2*sigma_L**2))*math.exp(-phi**2/(2*sigma_T**2))

"""" Chooses the model that best represents the transmission efficiency at dT
   Arguments:
      dT: a decimal representing the energy level of interest
    Returns:  
      A number indicating which model fits best:
      - 1 means gaussian
      - 2 means mod_gaussian
      - 3 means mod gaussian fixed sigma L
      - 0 means dT is between -0.15 (mod gauss) and -0.1 (gauss)
      - 4 means dT between 0(gauss) and 0.05(gauss fixed sigma L)""" 
def choose_model(dT):
  if dT <=-0.15:
    return 2
  elif -0.1<= dT<=0:
    return 1
  elif dT>= 0.05:
    return 3
  elif -0.15<dT<-0.1:
    return 0
  else:
    return 4
    
#these functions are only set up for discrete dT values
"""" Obtains the parameters for a dT that is best represented by a
  gaussian model
   Arguments:
      dT: a decimal representing the energy level of interest
    Returns:  
      Parameters to be used in the gaussian model
       when calculating the transmission efficiency for dT:
       sigma_L, sigma_R, sigma_B, sigma_T""" 
def gauss_param(dT):
  ##from Jay models powerpoint
  if dT == -0.1:
    sigma_L = 1.87308
    sigma_R = 9.5499
    sigma_B = 2.76533
    sigma_T = 1.42978
  elif dT == -0.05:
    sigma_L = 2.63716
    sigma_R = 32.9357
    sigma_B = 2.8416
    sigma_T = 1.83227
  else:
    sigma_L = 5.3558
    sigma_R = 3.56623
    sigma_B = 2.15078
    sigma_T = 1.86417
  return sigma_L, sigma_R, sigma_B, sigma_T

"""" Obtains the parameters for any dT
   Arguments:
      dT: a decimal representing the energy level of interest
    Returns:  
      Parameters to be used in the modified gaussian model
       when calculating the transmission efficiency for dT:
       sigma_L, sigma_R, sigma_B, sigma_T, theta0, N""" 
def mod_gauss_param(dT):
  ##from Jay models powerpoint
  if dT == -0.2:
    sigma_L = 0.263606
    sigma_R = 0.10956
    sigma_B = 49.8503
    sigma_T = 0.524778
    theta0 = 1.19031
    N = 0.999998
  elif dT == -0.15:
    sigma_L = 2.33371
    sigma_R = 1.48944
    sigma_B = 3.88653
    sigma_T = 0.804163
    theta0 = 1.43261
    N = 0.797831
  elif dT == -0.1:
    sigma_L = 3.38933
    sigma_R = 8.56057
    sigma_B = 3.71699
    sigma_T = 1.58143
    theta0 = 1.37556
    N = 0.946582
  elif dT == -0.05:
    sigma_L = 4.23033
    sigma_R = 13.5232
    sigma_B = 3.13257
    sigma_T = 1.92483
    theta0 = 1.15403
    N = 0.986219
  elif dT == 0:
    sigma_L = 32.079
    sigma_R = 21.7662
    sigma_B = 2.48789
    sigma_T = 2.09244
    theta0 = -8.15265
    N = 0.94805
  elif dT == 0.05:
    sigma_L = 8.30103
    sigma_R = 5.65914
    sigma_B = 2.11787
    sigma_T = 2.29035
    theta0 = -3.61506
    N = 0.998989
  elif dT == 0.1:
    sigma_L = 7.90543
    sigma_R = 3.93158
    sigma_B = 1.84238
    sigma_T = 2.78002
    theta0 = -3.24628
    N = 1
  elif dT == 0.15:
    sigma_L = 8.58356
    sigma_R = 3.07716
    sigma_B = 1.48695
    sigma_T = 2.4876
    theta0 = -3.0238
    N = 0.999993
  else:
    sigma_L = 10.9377
    sigma_R = 2.41959
    sigma_B = 1.19343
    sigma_T = 2.69817
    theta0 = -3.01734
    N = 1
  return sigma_L, sigma_R, sigma_B, sigma_T, theta0, N

"""" Obtains the parameters for a dT that is best represented by a
 modified gaussian fixed sigma L model
   Arguments:
      dT: a decimal representing the energy level of interest
    Returns:  
      Parameters to be used in the modified gaussian fixed sigma L model
       when calculating the transmission efficiency for dT:
       sigma_R, sigma_B, sigma_T, theta0""" 
def mod_gauss_fsl_param(dT):
  ##from Jay models powerpoint
  if dT == 0:
    sigma_R = 0.101991
    sigma_B = 2.22355
    sigma_T = 1.91523
    theta0 = 2.58663
  elif dT == 0.05:
    sigma_R = 5.72369
    sigma_B = 2.12203
    sigma_T = 2.29458
    theta0 = -3.68096
  elif dT == 0.1:
    sigma_R = 3.92879
    sigma_B = 1.84293
    sigma_T = 2.7908
    theta0 = -3.24569
  elif dT == 0.15:
    sigma_R = 3.0754
    sigma_B = 1.48632
    sigma_T = 2.48725
    theta0 = -3.02201
  else:
    sigma_R = 2.41954
    sigma_B = 1.19299
    sigma_T = 2.69763
    theta0 = -3.01753
  return sigma_R, sigma_B, sigma_T, theta0

"""" Obtains the parameters for a dT that is best represented by a specific model
   Arguments:
      model: a number representing which model is to be used
      dT: a decimal representing the energy level of interest
    Returns:  
      Parameters to be used in the model when calculating the transmission efficiency for dT
    Requirements: 
      model must be 1, 2 or 3""" 
def get_parameters(model,dT):
  if model ==1:
    return gauss_param(dT)
  elif model == 2:
    return mod_gauss_param(dT)
  else:
    return mod_gauss_fsl_param(dT)
  
"""" Interpolates the transmission efficiency of a dT between two others
   Arguments:
      dT: a decimal representing the energy level of interest
      lower_dT: the lower x-bound for the interpolation
      upper_dT: the upper x-bound for the interpolation
      lower_e: the effiency that corresponds to lower_dT
      upper_e: the effiency that corresponds to upper_dT
    Returns:  
      Transmission efficiency of dT""" 
def interpolate(dT, lower_dT, upper_dT, lower_e, upper_e):
  slope = (upper_e-lower_e)/(upper_dT-lower_dT)
  return lower_e+slope*(dT- lower_dT)

"""" Calculates the transmission efficiency of a triplet of which the energy level 
is best fitted by a gaussian (-0.1<dT<-0.05)
   Arguments:
      dT: a decimal representing the energy level
      theta: an angle in degrees representing the horizontal angle
      phi: angle in degrees representing the vertical angle
    Returns:  
      Transmission efficiency of this particular""" 
def use_gaussian(dT, theta,phi):
  model = 1
  if -0.1<=dT<-0.05:
      dT_lower = -0.1
      dT_upper = -0.05

      sigma_L, sigma_R, sigma_B, sigma_T = get_parameters(model, dT_lower)
      result_lower = gaussian(theta, phi, sigma_L, sigma_R, sigma_B, sigma_T)

      sigma_L, sigma_R, sigma_B, sigma_T = get_parameters(model, dT_upper)
      result_upper = gaussian(theta, phi, sigma_L, sigma_R, sigma_B, sigma_T)

      return interpolate(dT, dT_lower, dT_upper, result_lower, result_upper)
  elif -0.05<=dT<=0:
      dT_lower = -0.05
      dT_upper =0

      sigma_L, sigma_R, sigma_B, sigma_T = get_parameters(model, dT_lower)
      result_lower = gaussian(theta, phi, sigma_L, sigma_R, sigma_B, sigma_T)

      sigma_L, sigma_R, sigma_B, sigma_T = get_parameters(model, dT_upper)
      result_upper = gaussian(theta, phi, sigma_L, sigma_R, sigma_B, sigma_T)

      return interpolate(dT, dT_lower, dT_upper, result_lower, result_upper)

"""" Calculates the transmission efficiency of a triplet of any energy level using a modified gaussian
   Arguments:
      dT: a decimal representing the energy level
      theta: an angle in degrees representing the horizontal angle
      phi: angle in degrees representing the vertical angle
    Returns:  
      Transmission efficiency of this particular""" 
def use_mod_gaussian(dT, theta, phi):
  if dT == -0.2 or dT == -0.15 or dT == -0.1 or dT == -0.05 or dT == 0 or dT == 0.05 or dT == 0.1 or dT == 0.15 or dT == 0.2:
      sigma_L, sigma_R, sigma_B, sigma_T, theta0, N = get_parameters(2, dT)
      return mod_gaussian(theta, phi, sigma_L, sigma_R, sigma_B, sigma_T, theta0, N)
  else:
    if -0.2<dT<-0.15:
      dt_lower = -0.2
      dt_upper = -0.15

    elif -0.15<dT<-0.1:
      dt_lower = -0.15
      dt_upper = -0.1

    elif -0.1<dT<-0.05:
      dt_lower = -0.1
      dt_upper = -0.05

    elif -0.05<dT<0:
      dt_lower = -0.05
      dt_upper = 0

    elif 0<dT<0.05:
      dt_lower = 0
      dt_upper = 0.05

    elif 0.05<dT<0.1:
      dt_lower = 0.05
      dt_upper = 0.1

    elif 0.1<dT<0.15:
      dt_lower = 0.1
      dt_upper = 0.15

    else:
      dt_lower = 0.15
      dt_upper = 0.2

    sigma_L_l, sigma_R_l, sigma_B_l, sigma_T_l, theta0_l, N_l = get_parameters(2,dt_lower)
    result_lower = mod_gaussian(theta, phi, sigma_L_l, sigma_R_l, sigma_B_l, sigma_T_l, theta0_l, N_l)

    sigma_L, sigma_R, sigma_B, sigma_T, theta0, N = get_parameters(2, dt_upper)
    result_upper = mod_gaussian(theta, phi, sigma_L, sigma_R, sigma_B, sigma_T, theta0, N)
    return interpolate(dT, dt_lower, dt_upper, result_lower, result_upper)      

"""" Calculates the transmission efficiency of a triplet of which the energy level 
is best fitted by a modified gaussian fixed sigma L model (dT>0)
   Arguments:
      dT: a decimal representing the energy level
      theta: an angle in degrees representing the horizontal angle
      phi: angle in degrees representing the vertical angle
    Returns:  
      Transmission efficiency of this particular""" 
def use_mod_gauss_fsl(dT, theta,phi):
  model = 3
  if dT == 0.05 or dT == 0.1 or dT ==0.15 or dT==0.2:
    sigma_R, sigma_B, sigma_T, theta0 = get_parameters(model, dT)
    return mod_gaussian_fixed_sigmaL(theta, phi, sigma_R, sigma_B, sigma_T, theta0)
  else:
    if 0.05<dT<0.1:
      dt_lower = 0.05
      dt_upper = 0.1
      
    elif 0.1<dT<0.15:
      dt_lower = 0.1
      dt_upper = 0.15
      
    else:
      dt_lower = 0.15
      dt_upper = 0.2
    
    sigma_R, sigma_B, sigma_T, theta0 = get_parameters(model, dt_lower)
    result_lower = mod_gaussian_fixed_sigmaL(theta, phi, sigma_R, sigma_B, sigma_T, theta0)

    sigma_R, sigma_B, sigma_T, theta0 = get_parameters(model, dt_upper)
    result_upper = mod_gaussian_fixed_sigmaL(theta, phi, sigma_R, sigma_B, sigma_T, theta0)
    return interpolate(dT, dt_lower, dt_upper, result_lower, result_upper)

"""" Calculates the transmission efficiency of a triplet of which the energy level 
is between modified gaussian and gaussian models (-0.15<dT<-0.1)
   Arguments:
      dT: a decimal representing the energy level
      theta: an angle in degrees representing the horizontal angle
      phi: angle in degrees representing the vertical angle
    Returns:  
      Transmission efficiency of this particular""" 
def between_mod_and_gauss(dT,theta,phi):
  sigma_L, sigma_R, sigma_B, sigma_T, theta0, N = get_parameters(2, -0.15)
  result_mod_gauss = mod_gaussian(theta, phi, sigma_L, sigma_R, sigma_B, sigma_T, theta0, N)

  sigma_L, sigma_R, sigma_B, sigma_T = get_parameters(1, -0.1)
  result_gauss = gaussian(theta, phi, sigma_L, sigma_R, sigma_B, sigma_T)

  return interpolate(dT, -0.15, -0.1, result_mod_gauss, result_gauss)

"""" Calculates the transmission efficiency of a triplet of which the energy level 
is between gaussian and modified gaussian fixed sigma L models (-0.05<dT<0)
   Arguments:
      dT: a decimal representing the energy level
      theta: an angle in degrees representing the horizontal angle
      phi: angle in degrees representing the vertical angle
    Returns:  
      Transmission efficiency of this particular"""
def between_gauss_and_mgfsl(dT, theta, phi):
  sigma_L, sigma_R, sigma_B, sigma_T = get_parameters(1, 0)
  result_gauss = gaussian(theta, phi, sigma_L, sigma_R, sigma_B, sigma_T)

  sigma_R, sigma_B, sigma_T, theta0 = get_parameters(3, 0.05)
  result_mgauss_fsl = mod_gaussian_fixed_sigmaL(theta, phi, sigma_R, sigma_B, sigma_T, theta0)

  return interpolate(dT, 0, 0.05, result_gauss, result_mgauss_fsl)

"""" Finds the model (or models) to use and returns the result of this function
   (the transmission efficiency for the given triplet)
   Arguments:
      dT: a decimal representing the energy level
      theta: an angle in degrees representing the horizontal angle
      phi: angle in degrees representing the vertical angle
    Returns:  
      Transmission efficiency of this particular"""
def pick_model_and_calc(dT,theta,phi):
  alpha = np.arccos(math.cos(theta*pi/180)*math.cos(phi*pi/180)+math.sin(theta*pi/180)*math.sin(phi*pi/180))
  if alpha*180/pi > 4.2 or alpha*180/pi < -4.2:
    return 0
  model = choose_model(dT)
  if model == 1:
    return use_gaussian(dT, theta, phi)
  elif model == 2:
    return use_mod_gaussian(dT, theta, phi)
  elif model == 3:
    return use_mod_gauss_fsl(dT, theta, phi)
  elif model == 0:
    return between_mod_and_gauss(dT,theta,phi)
  else:
    return between_gauss_and_mgfsl(dT,theta,phi)

def pick_model_and_calc_trunc(dT,theta,phi):
  if theta>3 or theta<-3 or phi>3 or phi<-3:
    return 0
  else:
    return pick_model_and_calc(dT,theta,phi)
def calc_rel_measurement_error(dT):
  if dT<uf.iloc[0,0]:
    return uf["Full aperture relative uncertainty"][0]
  for i in uf.index:
    lower_dT = uf["dT"][i]
    lower_uncertainty = uf["Full aperture relative uncertainty"][i]
    if uf.iloc[-1,0] == uf.iloc[i,0]:
      return lower_uncertainty
    elif dT == lower_dT:
      return lower_uncertainty
    upper_dT = uf["dT"][i+1]
    upper_uncertainty = uf["Full aperture relative uncertainty"][i+1]
    if lower_dT<dT<upper_dT:
      slope = (upper_uncertainty-lower_uncertainty)/(upper_dT-lower_dT)
      return lower_uncertainty+slope*(dT- lower_dT)

extrapolate_e = []
model_errors = []
measurement_errors = []
trunc_e = []
for i in df.index:
  dT = df["dt"][i]
  theta = df["theta"][i]
  phi = df["phi"][i]
  result = pick_model_and_calc(dT,theta,phi)
  extrapolate_e.append(result)
  if dT > -0.15 and result != 0:
    result_2 = use_mod_gaussian(dT, theta, phi)
    rel_error = abs(result_2-result)/result
    model_errors.append(rel_error)
  measurement_errors.append(calc_rel_measurement_error(dT))
  trunc_e.append(pick_model_and_calc_trunc(dT,theta,phi))

model_rel_error = np.mean(model_errors)
measurement_rel_error = np.mean(measurement_errors)
efficiency = np.mean(extrapolate_e)
trunc_efficiency = np.mean(trunc_e)
avg_e = (efficiency+trunc_efficiency)/2
trunc_rel_error = (efficiency-trunc_efficiency)/6/efficiency
total_e_error = (model_rel_error**2+measurement_rel_error**2+trunc_rel_error**2)**(1/2)
print('Extrapolated transmission efficiency is:',efficiency)
print('Truncated transmission efficiency is:', trunc_efficiency)
print('Average efficiency is:',avg_e)
print('Relative truncation error is:', trunc_rel_error)
print('Relative model error is:', model_rel_error)
print('Relative measurement error is:', measurement_rel_error)
print('total error:',total_e_error)