From c18d5c08c57f3549ed2db611f5aef253cb63a4bb Mon Sep 17 00:00:00 2001 From: Scott Prahl Date: Wed, 13 Mar 2024 09:56:20 -0700 Subject: [PATCH] improve reading .rxt files --- iadpython/ad.py | 12 +- iadpython/iad.py | 44 +++---- iadpython/port.py | 6 +- iadpython/rxt.py | 145 +++++++++++++++++------ iadpython/sphere.py | 283 ++++++++++++++++++++++++++++---------------- 5 files changed, 324 insertions(+), 166 deletions(-) diff --git a/iadpython/ad.py b/iadpython/ad.py index 65a01ba..3932e3c 100644 --- a/iadpython/ad.py +++ b/iadpython/ad.py @@ -26,8 +26,16 @@ def stringify(form, x): if x is None: - return 'None' - return form % x + s = 'None' + elif np.isscalar(x): + s = form % x + else: + mn = min(x) + mx = max(x) + s = form % mn + s += ' to ' + s += form % mx + return s class Sample(): """Container class for details of a sample. diff --git a/iadpython/iad.py b/iadpython/iad.py index 1b4276b..0f514c8 100644 --- a/iadpython/iad.py +++ b/iadpython/iad.py @@ -32,6 +32,19 @@ import scipy.optimize import iadpython as iad +def stringify(form, x): + if x is None: + s = 'None' + elif np.isscalar(x): + s = form % x + else: + mn = min(x) + mx = max(x) + s = form % mn + s += ' to ' + s += form % mx + return s + class Experiment(): """Container class for details of an experiment.""" @@ -93,8 +106,14 @@ def __str__(self): s = "---------------- Sample ---------------\n" s += self.sample.__str__() s += "\n--------------- Spheres ---------------\n" - if self.num_spheres == 0: + if not np.isscalar(self.num_spheres): + s += "number of spheres range (%s)\n" % stringify("%d",self.num_spheres) + elif self.num_spheres == 0: s += "No spheres used.\n" + elif self.num_spheres == 1: + s += "A single integrating sphere was used.\n" + elif self.num_spheres == 2: + s += "Double integrating spheres were used.\n" if self.r_sphere is not None: s += "Reflectance Sphere--------\n" s += self.r_sphere.__str__() @@ -105,28 +124,11 @@ def __str__(self): if self.include_measurements: s += "\n------------- Measurements ------------\n" s += " Reflection = " - if self.m_r is None: - s += "Missing\n" - elif np.isscalar(self.m_r): - s += "%.5f\n" % self.m_r - else: - s += "%s\n" % self.m_r.__str__() - + s += stringify("%.5f", self.m_r) + "\n" s += " Transmission = " - if self.m_t is None: - s += "Missing\n" - elif np.isscalar(self.m_r): - s += "%.5f\n" % self.m_t - else: - s += "%s\n" % self.m_t.__str__() - + s += stringify("%.5f", self.m_t) + "\n" s += " Unscattered Transmission = " - if self.m_u is None: - s += "Missing\n" - elif np.isscalar(self.m_r): - s += "%.5f\n" % self.m_u - else: - s += "%s\n" % self.m_u.__str__() + s += stringify("%.5f", self.m_u) + "\n" return s def check_measurements(self): diff --git a/iadpython/port.py b/iadpython/port.py index f6c6016..65743c0 100644 --- a/iadpython/port.py +++ b/iadpython/port.py @@ -187,9 +187,9 @@ def hit(self): r2 = (self.x - self.sphere.x)**2 r2 += (self.y - self.sphere.y)**2 r2 += (self.z - self.sphere.z)**2 - print('cap center (%7.2f, %7.2f, %7.2f)' % (self.x,self.y,self.z)) - print('pt on sph (%7.2f, %7.2f, %7.2f)' % (self.sphere.x,self.sphere.y,self.sphere.z)) - print("cap distance %7.2f %7.2f"% (r2, self.chord2)) +# print('cap center (%7.2f, %7.2f, %7.2f)' % (self.x,self.y,self.z)) +# print('pt on sph (%7.2f, %7.2f, %7.2f)' % (self.sphere.x,self.sphere.y,self.sphere.z)) +# print("cap distance %7.2f %7.2f"% (r2, self.chord2)) return r2 < self.chord2 def set_center(self, x, y, z): diff --git a/iadpython/rxt.py b/iadpython/rxt.py index 6125a1c..0363be1 100644 --- a/iadpython/rxt.py +++ b/iadpython/rxt.py @@ -15,7 +15,7 @@ >>> import iadpython >>> >>> filename = 'ink.rxt' - >>> exp = iadpython.read_iad_input(filename) + >>> exp = iadpython.read_rxt(filename) >>> if exp.lambda0 is None: >>> plt.plot(exp.m_r) >>> else: @@ -25,6 +25,7 @@ >>> plt.show() """ +import os import re import numpy as np import iadpython @@ -35,10 +36,13 @@ def read_and_remove_notation(filename): """Read file and remove all whitespace and comments.""" s = '' + + if not os.path.exists(filename): + raise ValueError('input file "%s" must end in ".rxt"' % filename) + with open(filename, encoding="utf-8") as f: for line in f: - line = re.sub(r'#.*', '', line) - line = re.sub(r'\s', ' ', line) + line = re.sub(r'\s*#.*', '', line) line = re.sub(r',', ' ', line) s += line @@ -47,10 +51,87 @@ def read_and_remove_notation(filename): s = re.sub(r'IAD1', '', s) s = re.sub(r'\s+', ' ', s) - s = re.sub(r'^\s+', '', s) - s = re.sub(r'\s+$', '', s) + s = s.rstrip() + s = s.lstrip() return s +def fill_in_data(exp, data, column_letters_str): + if column_letters_str == '': + columns = int(exp.num_measures) + if data[0]>1: + columns+=1 + else : + columns = len(column_letters_str) + + data_in_columns = data.reshape(-1, columns) + exp.lambda0 = None + if column_letters_str == '': + col = 0 + if data[0] > 1: + exp.lambda0 = data_in_columns[:,0] + col = 1 + exp.m_r = data_in_columns[:,col] + if exp.num_measures >=2: + exp.m_t = data_in_columns[:,col+1] + if exp.num_measures >=3: + exp.m_u = data_in_columns[:,col+2] + if exp.num_measures >=4: + exp.r_sphere.r_wall = data_in_columns[:,col+3] + if exp.num_measures >=5: + exp.t_sphere.r_wall = data_in_columns[:,col+4] + if exp.num_measures >=6: + exp.r_sphere.r_std = data_in_columns[:,col+5] + if exp.num_measures >=7: + exp.t_sphere.r_std = data_in_columns[:,col+6] + if exp.num_measures >7: + raise ValueError('unimplemented') + return + + for col, letter in enumerate(column_letters_str): + if letter=='a': + exp.default_a = data_in_columns[:,col] + elif letter=='b': + exp.default_b = data_in_columns[:,col] + elif letter=='B': + exp.d_beam = data_in_columns[:,col] + elif letter=='c': + exp.fraction_of_rc_in_mr = data_in_columns[:,col] + elif letter=='C': + exp.fraction_of_tc_in_mt = data_in_columns[:,col] + elif letter=='d': + exp.sample.d = data_in_columns[:,col] + elif letter=='D': + exp.sample.d_above = data_in_columns[:,col] + exp.sample.d_below = data_in_columns[:,col] + elif letter=='e': + exp.error = data_in_columns[:,col] + elif letter=='g': + exp.default_g = data_in_columns[:,col] + elif letter=='t': + exp.m_t = data_in_columns[:,col] + elif letter=='L': + exp.lambda0 = data_in_columns[:,col] + elif letter=='n': + exp.n = data_in_columns[:,col] + elif letter=='N': + exp.n_above = data_in_columns[:,col] + exp.n_below = data_in_columns[:,col] + elif letter=='r': + exp.m_r = data_in_columns[:,col] + elif letter=='R': + exp.r_sphere.r_std = data_in_columns[:,col] + elif letter=='S': + exp.num_spheres = data_in_columns[:,col] + elif letter=='T': + exp.t_sphere.r_rstd = data_in_columns[:,col] + elif letter=='u': + exp.m_u = data_in_columns[:,col] + elif letter=='w': + exp.r_sphere.r_wall = data_in_columns[:,col] + elif letter=='W': + exp.t_sphere.r_wall = data_in_columns[:,col] + else: + raise ValueError('unimplemented column type "%s"' % letter) def read_rxt(filename): """Read an IAD input file in .rxt format. @@ -62,7 +143,14 @@ def read_rxt(filename): Experiment object """ s = read_and_remove_notation(filename) - x = np.array([float(value) for value in s.split(' ')]) + x = s.split(' ') + + # Remove single-letter entries and save them + column_letters = [item for item in x if len(item) == 1 and item.isalpha()] + x = [item for item in x if not (len(item) == 1 and item.isalpha())] + column_letters_str = ''.join(column_letters) + + x = np.array([float(item) for item in x]) sample = iadpython.Sample(a=None, b=None, g=None) sample.n = x[0] @@ -93,35 +181,18 @@ def read_rxt(filename): if exp.num_spheres > 0: exp.t_sphere = iadpython.Sphere(x[12], x[13], x[14], x[15], 0, x[16]) - exp.num_measures = x[17] - - exp.lambda0 = np.zeros(0) - if exp.num_measures >= 1: - exp.m_r = np.zeros(0) - - if exp.num_measures >= 2: - exp.m_t = np.zeros(0) - - if exp.num_measures >= 3: - exp.m_u = np.zeros(0) - - count_per_line = 1 - for i in range(18, len(x)): - - if x[i] > 1: - exp.lambda0 = np.append(exp.lambda0, x[i]) - continue - if count_per_line == 1: - exp.m_r = np.append(exp.m_r, x[i]) - elif count_per_line == 2: - exp.m_t = np.append(exp.m_t, x[i]) - elif count_per_line == 3: - exp.m_u = np.append(exp.m_u, x[i]) - if count_per_line >= exp.num_measures: - count_per_line = 1 - else: - count_per_line += 1 - - if len(exp.lambda0) == 0: - exp.lambda0 = None + exp.num_measures = 0 + if column_letters_str == '': + exp.num_measures = x[17] + data = x[18:] + else: + if 'r' in column_letters_str: + exp.num_measures += 1 + if 't' in column_letters_str: + exp.num_measures += 1 + if 'u' in column_letters_str: + exp.num_measures += 1 + data = x[17:] + + fill_in_data(exp, data, column_letters_str) return exp diff --git a/iadpython/sphere.py b/iadpython/sphere.py index e4a9e6e..38a18bd 100644 --- a/iadpython/sphere.py +++ b/iadpython/sphere.py @@ -32,8 +32,7 @@ Methods: cap_area: actual area of a spherical cap for a port relative_cap_area: relative area of spherical cap to the sphere's area. - gain: sphere gain relative to a black sphere - multiplier: multiplier for wall power due to sphere + gain: sphere gain relative to isotropic diffuse source in center of black sphere Example usage: >>> import iadpython @@ -48,30 +47,29 @@ import random import numpy as np +from enum import Enum import iadpython -def uniform_disk(): - """ - Generate a point uniformly distributed on a unit disk. - - This function generates and returns a point (x, y) that is uniformly distributed - over the area of a unit disk centered at the origin (0, 0). The unit disk is - defined as the set of all points (x, y) such that x^2 + y^2 <= 1. The method - uses rejection sampling, where random points are generated within the square - that bounds the unit disk, and only points that fall within the disk are accepted. - - Returns: - tuple of (float, float, float): A point (x, y) where x and y are the - coordinates of the point uniformly distributed within the unit disk, - and s is the square of the distance from the origin to this point. - """ - s = 2 - while s > 1: - x = 2 * random.random() - 1 - y = 2 * random.random() - 1 - s = x*x + y*y - return x, y, s - +class PortType(Enum): + """Possible sphere wall locations.""" + EMPTY = 0 + WALL = 1 + SAMPLE = 2 + DETECTOR = 3 + FIRST = 4 + +def stringify(form, x): + if x is None: + s = 'None' + elif np.isscalar(x): + s = form % x + else: + mn = min(x) + mx = max(x) + s = form % mn + s += ' to ' + s += form % mx + return s class Sphere(): """Container class for an integrating sphere. @@ -123,74 +121,120 @@ def __str__(self): """Return basic details as a string for printing.""" s = "" s += "Sphere\n" - s += " diameter = %7.2f mm\n" % self.d - s += " radius = %7.2f mm\n" % (self.d/2) + s += " diameter = %s mm\n" % stringify("%7.2f", self.d) + s += " radius = %s mm\n" % stringify("%7.2f", self.d/2) s += " relative area = %7.4f\n" % self.a_wall - s += " uru walls = %7.1f%%\n" % (self.r_wall * 100) - s += " uru standard = %7.1f%%\n" % (self.r_std * 100) + s += " uru walls = %s\n" % stringify("%7.1f%%", self.r_wall * 100) + s += " uru standard = %s\n" % stringify("%7.1f%%", self.r_std * 100) s += " baffle = %s\n" % self.baffle s += "Sample Port\n" + str(self.sample) s += "Detector Port\n" + str(self.detector) s += "Empty Port\n" + str(self.empty) - s += "Multipliers\n" - s += " nothing = %7.3f\n" % self.multiplier(0, 0) - s += " standard = %7.3f\n" % self.multiplier(self.r_std, self.r_std) - s += "Monte Carlo\n" - s += " location = (%6.1f, %6.1f, %6.1f) mm\n" % (self.x, self.y, self.z) - s += " weight = %7.3f\n" % self.weight + s += "Gain range\n" + s += " nothing = %s\n" % stringify("%7.3f", self.gain(0.0)) + s += " standard = %s\n" % stringify("%7.3f", self.gain(self.r_std)) + s += " 100%% = %s\n" % stringify("%7.3f", self.gain(1.0)) return s - def multiplier(self, UX1=None, URU=None): + def gain(self, sample_uru=None): """ - Determine detected light in a sphere. - - The idea here is that UX1 is the diffuse light entering the sphere. This - might be diffuse light reflected by the sample or diffuse light transmitted - through the sample. - - We assume that the sample is characterized by the reflection and tranmission - for collimated normally incident light (UR1 and UT1) and by the reflection - and transmission for diffuse incident light (URU and UTU). - - There are four common cases - - 1. For a reflection experiment, `UX1=UR1` and `URU=URU` - - 2. For a transmission experiment, `UX1=UT1` and `URU=URU` - - 3. If light hits the sphere wall first, then `UX1=r_wall` - - 4. If light is enters the sphere completely diffuse then `UX1=1` - - This matches LabSphere, "Technical Guide: integrating Sphere Theory - and application" using equation 14 for case 3 + Determine gain due to multiple bounces in the sphere. + + If UX1 is the power passing through the sample UT1 or reflected by the + sample UR1, then power falling on the detector will be + + P_detector = a_detector * gain * UX1 * P_0 - Args: - UX1: diffuse light entering the sphere - URU: sample reflectance for diffuse irradiance - Returns: - sphere multiplier + and power detected will be + + P_detected = (1-r_detector) * a_detector * gain * UX1 * P_0 """ - if UX1 is None: - UX1 = self.r_wall - - if URU is None: - URU = UX1 - - denom = 1 - denom -= self._a_wall * self.r_wall - denom -= self.sample.a * URU - denom -= self.detector.a * self.detector.uru + if sample_uru is not None: + original_uru = self.sample.uru + self.sample.uru = sample_uru + + tmp = self.detector.a * self.detector.uru + self.sample.a * self.sample.uru - if np.isscalar(denom): - if denom > 1e-8: - m = UX1 * self._a_wall / denom - else: - m = np.inf + if self.baffle: + denom = 1 - self.r_wall * (self.a_wall + (1 - self.empty.a) * tmp) else: - m = np.where(denom > 1e-8, UX1 * self._a_wall / denom, np.inf) - return m - + denom = 1 - self._a_wall * self.r_wall - tmp + + g = 1 / denom + + # restore value + if sample_uru is not None: + self.sample.uru = original_uru + return g + + def pdetector(self): + P = (1-self.empty.a) * self.r_wall + N = 1000 + pd = np.zeros(N) + pw = np.zeros(N) + ps = np.zeros(N) + + pd[0] = self.detector.a * P + ps[0] = self.sample.a * P + pw[0] = self.a_wall * P + + for j in range(N-1): + pd[j+1] = self.detector.a * self.r_wall * pw[j] + ps[j+1] = self.sample.a * self.r_wall * pw[j] + pw[j+1] = self.a_wall * self.r_wall * pw[j] + pw[j+1] += (1 - self.empty.a) * self.detector.uru * pd[j] + pw[j+1] += (1 - self.empty.a) * self.sample.uru * ps[j] + + sumw = np.cumsum(pw) + sumw -= sumw[1] + sumd = np.cumsum(pd) + print(' k P_d^k P_w^k sum(P_d^k) sum(P_w^k)') + for j in range(10): + print("%3d %9.5f %9.5f %9.5f %9.5f" % (j+1, pd[j], pw[j], sumd[j], sumw[j])) + print("%3d %9.5f %9.5f %9.5f %9.5f" % (N-1, pd[N-1], pw[N-1], sumd[N-1], sumw[N-1])) + print() + + pd[0] = self.detector.a * P + ps[0] = self.sample.a * P + pw[0] = self.a_wall * P + + pw[1] = self.a_wall * self.r_wall * pw[0] + pw[1] += (1 - self.empty.a) * self.detector.uru * pd[0] + pw[1] += (1 - self.empty.a) * self.sample.uru * ps[0] + pd[1] = self.detector.a * self.r_wall * pw[0] + + beta = (1 - self.empty.a) + beta *= self.detector.a * self.detector.uru + self.sample.a * self.sample.uru + for j in range(1,N-1): + pw[j+1] = self.r_wall * (self.a_wall * pw[j] + beta * pw[j-1]) + + sumw = np.cumsum(pw) + sumw -= sumw[1] + sumd = np.cumsum(pd) + print(' k P_d^k P_w^k sum(P_d^k) sum(P_w^k)') + for j in range(10): + print("%3d %9.5f %9.5f %9.5f %9.5f" % (j+1, pd[j], pw[j], sumd[j], sumw[j])) + print("%3d %9.5f %9.5f %9.5f %9.5f" % (N-1, pd[N-1], pw[N-1], sumd[N-1], sumw[N-1])) + print() + + beta = (1 - self.empty.a) + beta *= self.detector.a * self.detector.uru + self.sample.a * self.sample.uru + numer = self.a_wall**3*self.r_wall+ beta*(2*self.a_wall + self.a_wall**2*self.r_wall + beta) + denom = 1-self.r_wall *(self.a_wall + beta) + sum3 = self.r_wall * numer/denom * P + + pdx = self.detector.a * P + pdx += self.detector.a * self.r_wall * self.a_wall * P + pdx += self.detector.a * self.r_wall * (self.a_wall**2 * self.r_wall + beta) * P + pdx += self.detector.a * self.r_wall * sum3 + + print("%9.5f" % pdx) + + beta = (1 - self.empty.a) + beta *= self.detector.a * self.detector.uru + self.sample.a * self.sample.uru + pdx = self.detector.a * P/(1-self.r_wall*(self.a_wall+beta)) + print("%9.5f" % pdx) + @property def d(self): """Getter property for sphere diameter.""" @@ -234,7 +278,7 @@ def r_std(self, value): else: assert 0 <= value.all() <= 1, "Reflectivity of standard must be between 0 and 1" self._r_std = value - self.gain_std = self.multiplier(self.r_std, self.r_std) + self.gain_std = self.gain(self.r_std) @property def r_wall(self): @@ -254,18 +298,21 @@ def uniform(self): """ Generate a point uniformly distributed on the surface of a sphere. - This function returns a point that is uniformly distributed over the surface - of a sphere with a specified radius. It utilizes the method described in - "Choosing a Point from the Surface of a Sphere" by George Marsaglia (1972), - which is an efficient algorithm for generating such points using a transformation - from points uniformly distributed on a unit disk. - + Using Gaussian distribution for all three coordinates ensures a + uniform distribution on the surface of the sphere. See + + https://math.stackexchange.com/questions/1585975 + Returns: - A numpy array of the coordinates (x, y, z) of a random point on the sphere's surface. + (x, y, z) for a random point on the sphere's surface. """ - ux, uy, s = uniform_disk() - root = 2 * np.sqrt(1-s) - return np.array([ux * root, uy * root, 1 - 2 * s]) * self.d/2 + while True: + x = random.gauss() + y = random.gauss() + z = random.gauss() + r = np.sqrt(x**2 + y**2 + z**2) + if r > 0: + return np.array([x, y, z]) * (self.d / 2) / r def do_one_photon(self): """Bounce photon inside sphere until it leaves.""" @@ -275,9 +322,12 @@ def do_one_photon(self): # assume photon launched form sample weight = 1 last_location = PortType.SAMPLE - +# R = self.d/2 while weight > 1e-4: + lastx = self.x + lasty = self.y + lastz = self.z self.x, self.y, self.z = self.uniform() if self.detector.hit(): @@ -287,7 +337,15 @@ def do_one_photon(self): continue # record detected light and update weight - transmitted = weight * (1-self.r_detector) + vx=self.x-lastx + vy=self.y-lasty + vz=self.z-lastz +# RR = np.sqrt(vx*vx+vy*vy+vz*vz) +# print("n = [%5.2f, %5.2f, %5.2f]" % (self.x/R, self.y/R, self.z/R)) +# print("v = [%5.2f, %5.2f, %5.2f]" % (vx/RR, vy/RR, vz/RR)) +# costheta = abs((vx * self.x) + (vy * self.y) + (vz * self.z))/R/RR +# print(costheta) + transmitted = weight * (1-self.detector.uru) detected += transmitted weight -= transmitted last_location = PortType.DETECTOR @@ -297,7 +355,7 @@ def do_one_photon(self): continue if last_location == PortType.DETECTOR and self.baffle: # detector --> sample prohibited continue - weight *= self.r_sample + weight *= self.sample.uru last_location = PortType.SAMPLE elif self.empty.hit(): @@ -313,18 +371,37 @@ def do_one_photon(self): return detected, bounces - def do_N_photons(self): + def do_N_photons(self, N): + + num_trials = 10 + total_detected = np.zeros(num_trials) + total_bounces = np.zeros(num_trials) + + N_per_trial = N // num_trials - total_detected = 0 - total_bounces = 0 + total = 0 + for j in range(num_trials): + for i in range(N_per_trial): + detected, bounces = self.do_one_photon() + # print("%d %8.3f %d" % (i,detected, bounces)) + total_detected[j] += detected + total_bounces[j] += bounces + total += 1 - for i in range(N): - detected, bounces = self.do_one_photon() - total_detected += detected - total_bounces += bounces + ave = np.mean(total_detected)/N_per_trial + std = np.std(total_detected)/N_per_trial + stderr = std / np.sqrt(num_trials) + print("average detected = %.3f ± %.3f" % (ave,stderr)) + ave /= self.detector.a * (1-self.detector.uru) + std /= self.detector.a * (1-self.detector.uru) + stderr = std / np.sqrt(num_trials) + print("average gain = %.3f ± %.3f" % (ave,stderr)) + print("calculated gain = %.3f" % self.gain()) - print("average detected = %.5f" % (total_detected/N)) - print("average bounces = %.5f" % (total_bounces/N)) + ave = np.mean(total_bounces)/N_per_trial + std = np.std(total_bounces)/N_per_trial + stderr = std / np.sqrt(num_trials) + print("average bounces = %.3f ± %.3f" % (ave,stderr)) def Gain_11(RS, TS, URU, tdiffuse):