import numpy as npy from typing import Dict import scipy.special import pandas as pd import PyPWA as pwa from dataclasses import dataclass import numexpr as ne ## PWA Simulation Tutorial (Template)## ## Defines a class with functions that are introduced to PyPWA through pwa.NestedFunction @dataclass class Resonances: resonance_index: int wave_index: int Cr: npy.ndarray Wave: npy.ndarray wave_data: npy.ndarray W0: npy.ndarray R0: npy.ndarray class NewAmplitude(pwa.NestedFunction): ## Necessary functions: "setup" and "calculate" def setup(self, kVars): self.__kVars = kVars self.__res = self.__load_resonance() self.__Vs = self.__make_Vs() self.__decay = self.__produce_specific_decay(self.__res.wave_data) self.__resonance = None def calculate(self, *args): intensity = npy.zeros(len(self.__kVars), "c16") rho = npy.zeros(len(self.__kVars), "c16") for res1 in range(self.__res.resonance_index): for res2 in range(self.__res.resonance_index): for wave1 in range(self.__res.wave_index): for wave2 in range(self.__res.wave_index): rho.real=0 rho.imag=0 Vs2 = self.__Vs[res1][wave1] * self.__Vs[res2][wave2].conjugate() decay2 = self.__decay[wave1] * self.__decay[wave2].conjugate() #Unpolarized Intensity # self.__kVars["pol"]= 0. if self.__res.wave_data["E"][wave1] == 1 and self.__res.wave_data["E"][wave2] == 1: rho.real = 1 - self.__kVars["pol"]*npy.cos(2*self.__kVars["alpha"]) if self.__res.wave_data["E"][wave1] == 1 and self.__res.wave_data["E"][wave2] == -1: rho.imag = -1 * self.__kVars["pol"]*npy.sin(2 * self.__kVars["alpha"]) if self.__res.wave_data["E"][wave1] == -1 and self.__res.wave_data["E"][wave2] == 1: rho.imag = 1 * self.__kVars["pol"]*npy.sin(2 * self.__kVars["alpha"]) if self.__res.wave_data["E"][wave1] == -1 and self.__res.wave_data["E"][wave2] == -1: rho.real = 1 + self.__kVars["pol"]*npy.cos(2*self.__kVars["alpha"]) intensity += 0.5 * (Vs2 * decay2) * rho # Kill's the intensity if there is still imaginary values left, and should display # those imaginary values if npy.sum(intensity.imag) > 1.e-6: print("Results are wrong! Imaginary values have been found!") self.error = intensity return intensity.real #-------------------------------------------------------------------- ## Auxiliary functions and definitions ## Define resonances and waves ************** def __load_resonance(self): res = pd.read_csv("resonanceNEW3.csv").to_records(index=False) elm = pd.read_csv("wavesNEW3.csv").to_records(index=False) res=pd.DataFrame(res) elm["E"] = elm["E"].astype(int) elm["ER"] = elm["ER"].astype(int) elm["L"] = elm["L"].astype(int) elm["M"] = elm["M"].astype(int) return self.__array_to_resonance(res, elm) def __make_Vs(self): Vs = npy.empty([self.__res.resonance_index, self.__res.wave_index, len(self.__kVars)], "c16") for index in range(self.__res.resonance_index): for w_index in range(self.__res.wave_index): Cr = self.__res.Cr[index] Wave = self.__res.Wave[index][w_index] W0 = self.__res.W0[index] R0 = self.__res.R0[index] mass = self.__kVars["mass"] frac = ne.evaluate("(sqrt(Cr * Wave)*W0*R0)/((W0**2 - mass**2)**2 + W0**2 * R0**2)") complex_part = ne.evaluate("(W0**2 - mass**2) + 1j*(W0 * R0)") Vs[index][w_index] = frac * complex_part return Vs def __array_to_resonance(self, array: npy.ndarray, elm: npy.ndarray): # res_index = len(array) # wave_index = len(array.dtype) - 3 r,c = array.shape res_index=r wave_index=c-3 Cr = array["Cr"].astype("f4") W0 = array["W0"].astype("f4") R0 = array["R0"].astype("f4") Wave = npy.empty([res_index, wave_index], "f4") for index in range(res_index): for w_index in range(wave_index): Wave[index][w_index] = array[f"Wave{w_index}"][index] self.__resonance = Resonances(res_index, wave_index, Cr, Wave, elm, W0, R0) return self.__resonance #Production Amplitudes def __make_Vs_old(self): Vs = npy.empty([self.__res.resonance_index, self.__res.wave_index, len(self.__kVars)], "c16") complex_part = npy.empty(len(self.__kVars), "c16") for index in range(self.__res.resonance_index): for w_index in range(self.__res.wave_index): numerator = npy.sqrt(self.__res.Cr[index] * self.__res.Wave[index][w_index]) numerator *= self.__res.W0[index] * self.__res.R0[index] denominator = (self.__res.W0[index]**2 - self.__kVars["mass"]**2)**2 denominator += self.__res.W0[index]**2 * self.__res.R0[index]**2 # denominator /= self.__norm complex_real = self.__res.W0[index]**2 - self.__kVars["mass"]**2 complex_imag = self.__res.W0[index] * self.__res.R0[index] complex_part.real = complex_real complex_part.imag = complex_imag Vs[index][w_index] = (numerator/denominator) * complex_part return Vs #Decay Amplitudes def __produce_specific_decay(self, wave_data): thetaC = self.__kVars["phi"].copy() phiC = self.__kVars["theta"].copy() thetaC[thetaC<0] = thetaC[thetaC<0]+2*npy.pi decay = npy.empty([len(wave_data), len(self.__kVars)], "c16") for index, wave in enumerate(wave_data): decay[index] = scipy.special.sph_harm(wave["M"], wave["L"], thetaC, phiC) decay[index] -= wave["ER"]*npy.power((-1),wave["M"])*scipy.special.sph_harm(-wave["M"], wave["L"], thetaC, phiC) return decay def __norm_int(self): for wave1 in range(self.__res.wave_index): for wave2 in range(self.__res.wave_index): decay = self.__decay[wave1] * self.__decay[wave2].conjugate() if self.__res.wave_data["E"][wave1] == 1 and self.__res.wave_data["E"][wave2] == 1: rho = 1-self.__kVars["pol"]*npy.cos(2*self.__kVars["alpha"]) if self.__res.wave_data["E"][wave1] == 1 and self.__res.wave_data["E"][wave2] == -1: rho = -1j * self.__kVars["pol"]*npy.sin(2 * self.__kVars["alpha"]) if self.__res.wave_data["E"][wave1] == -1 and self.__res.wave_data["E"][wave2] == 1: rho = 1j * self.__kVars["pol"] * npy.sin(2 * self.__kVars["alpha"]) if self.__res.wave_data["E"][wave1] == -1 and self.__res.wave_data["E"][wave2] == -1: rho = 1 + self.__kVars["pol"] * npy.cos(2 * self.__kVars["alpha"]) # intensity += .5 * Vs * decay * rho # no considering polarization normInt += decay normInt /= len(decay) return normInt @property def Vs(self): return self.__Vs.copy() @Vs.setter def Vs(self, variable): if len(variable) != len(self.__Vs): raise ValueError(f"Must be {len(self.__Vs)} length!") else: self.__Vs = variable @property def decay(self): return self.__decay.copy() @decay.setter def decay(self, new_decay): self.__decay = new_decay @property def resonance(self): return self.__res