Compute integration from goniophotometer measurements.

get_pgii_tdd.py

import argparse
import glob
import numpy as np
import scipy.spatial


def trapezoidal_area(xyz):
    """Calculate volume under a surface defined by irregularly spaced points
    using delaunay triangulation. "x,y,z" is a <numpoints x 3> shaped ndarray."""
    grid = scipy.spatial.Delaunay(xyz[:,:2])
    tri = xyz[grid.vertices]
    a = tri[:,0,:2] - tri[:,1,:2]
    b = tri[:,0,:2] - tri[:,2,:2]
    vol = np.cross(a, b) @ tri[:,:,2]
    return vol.sum() / 6.0


def get_data(dirpath, inthe, inphi):
    """Get the relevant files."""
    files = []
    for fpath in glob.glob(dirpath):
        with open(fpath) as rdr:
            row = rdr.readline()
            _intheta = 0
            _inphi = 0
            while row.startswith('#'):
                _key_val = row.split()
                _key = _key_val[0]
                if _key == '#intheta':
                    _intheta = float(_key_val[1])
                if _key == '#inphi':
                    _inphi = float(_key_val[1])
                row = rdr.readline()
            if _intheta == inthe and _inphi == inphi:
                files.append(fpath)
    if len(files) > 1:
        raise ValueError("More than one file has the requested incidence")
    if len(files) == 0:
        raise ValueError("Cannot find requested incident angle")
    return np.loadtxt(files[0])


def get_filtered_tdd(data, tin, pin, ha,
                     dsf=True, refl=False):
    xpos = np.sin(np.radians(data[:,0])) * np.cos(np.radians(data[:,1]))
    ypos = np.sin(np.radians(data[:,0])) * np.sin(np.radians(data[:,1]))
    zpos = np.cos(np.radians(data[:,0]))
    if refl:
        pin += 180
    else:
        tin += 180
    xin = np.sin(np.radians(tin)) * np.cos(np.radians(pin))
    yin = np.sin(np.radians(tin)) * np.sin(np.radians(pin))
    zin = np.cos(np.radians(tin))
    oac = np.cos(np.radians(ha))
    angles = xpos * xin + ypos * yin + zpos * zin
    if ha == 90:
        filtered = data
    else:
        filtered = data[np.where(angles > oac)]
    xyz = np.zeros(filtered.shape)
    xyz[:,0] = np.sin(np.radians(filtered[:,0])) * np.cos(np.radians(filtered[:,1]))
    xyz[:,1] = np.sin(np.radians(filtered[:,0])) * np.sin(np.radians(filtered[:,1]))
    if dsf:
        xyz[:,2] = filtered[:,2]/np.abs(np.cos(np.radians(filtered[:,0])))
    else:
        xyz[:,2] = filtered[:,2]
    return trapezoidal_area(xyz)

if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument('-f', required=True)
    parser.add_argument('-tin', type=float, required=True)
    parser.add_argument('-pin', type=float, required=True)
    parser.add_argument('-ha', type=float, required=True, help="Tdir-dir half angle")
    parser.add_argument('-r', action='store_true')
    parser.add_argument('-cos', action='store_false')
    args = parser.parse_args()
    data = get_data(args.f, args.tin, args.pin)
    integral = get_filtered_tdd(data, args.tin, args.pin, args.ha,
                                dsf=args.cos, refl=args.r)
    print(integral)