'''Calculate the integrated irradiance from one or more StellarNet irradiance files. The integral is calculated using the trapezoidal rule. ''' from __future__ import division import string, sys, getopt from math import log10 # The following weighting curve came from Appendix B of # http://www.hpa.org.uk/radiation/publications/documents_of_nrpb/abstracts/absd13-3.htm # # It is used to calculate the effective irradiance of UV radiation # (see paragraph 33, page 7 in the paper): # # "The basic exposure limit value (ELV) for both general public and # occupational exposure to UVR (ultraviolet radiation) incident # perpendicular to the skin or eye is an effective radiant exposure of # 30 J m^(-2) within any 8-hour period (ICNIRP, 1999). This assumes # that the exposure is delivered during any period of 8 successive # hours, even where such a period overlaps work shifts or calendar # days. Effective radiant exposure (J m^(-2)) is the product of # exposure duration (s) and effective irradiance (W m^(-2)). # Effective irradiance is the spectral irradiance Elambda at the eye # or skin surface, mathematically weighted with the hazard relative # spectral effectiveness factor s(lambda) and summed from 180 to 400 # nm, as follows: # # Eeff = SUM[Elambda*s(lambda)*delta_lambda] # # where # Eeff = effective irradiance in W m^(-2) # Elambda = spectral irradiance from measurements in W m^-2 nm^-1 # s(lambda) = relative spectral effectiveness factor # delta_lambda = bandwidth of the calculation or measurement in nm" # # The relative spectral effectiveness factor is given in the following # array: weighting_curve = ( (180, 0.012), (190, 0.019), (200, 0.030), (205, 0.051), (210, 0.075), (215, 0.095), (220, 0.120), (225, 0.150), (230, 0.190), (235, 0.240), (240, 0.300), (245, 0.360), (250, 0.430), (254, 0.500), (255, 0.520), (260, 0.650), (265, 0.810), (270, 1.000), (275, 0.960), (280, 0.880), (285, 0.770), (290, 0.640), (295, 0.540), (297, 0.460), (300, 0.300), (303, 0.120), (305, 0.060), (308, 0.026), (310, 0.015), (313, 0.006), (315, 0.003), (316, 0.0024), (317, 0.0020), (318, 0.0016), (319, 0.0012), (320, 0.0010), (322, 0.00067), (323, 0.00054), (325, 0.00050), (328, 0.00044), (330, 0.00041), (333, 0.00037), (335, 0.00034), (340, 0.00028), (345, 0.00024), (350, 0.00020), (355, 0.00016), (360, 0.00013), (365, 0.00011), (370, 0.000093), (375, 0.000077), (380, 0.000064), (385, 0.000053), (390, 0.000044), (395, 0.000036), (400, 0.000030), ) # I fitted straight lines to the above curve to allow analytical # evaluation. The lines are between the following points and # undefined outside of them. Note the y values need to be the log, as # the lines are for a semilog plot. weighting_points = ( (180, 0.012, 200, 0.03), (200, 0.030, 210, 0.075), (210, 0.075, 270, 1.000), (270, 1.000, 295, 0.600), (295, 0.600, 323, 5.4e-4), (323, 5.4e-4, 400, 3e-5) ) weighting = None # Will contain linear approximation constants calculate_weighted = 0 # Use weighting for irr if true def Weighting(lambda1): '''Return the weighting factor for wavelengths between 180 nm and 400 nm. Return 1e-99 if outside this range. ''' global weighting if weighting == None: weighting = [] for la1, w1, la2, w2 in weighting_points: slope = (log10(w2) - log10(w1))/(la2 - la1) intercept = log10(w1) - slope*la1 weighting.append((la1, la2, slope, intercept)) weighting = tuple(weighting) for la1, la2, slope, intercept in weighting: if la1 <= lambda1 <= la2: y = slope*lambda1 + intercept return 10**y return 1e-99 def TestWeighting(): '''Print out the ratio of the calculated vs actual. Around 300-320 nm, there are 20-30% deviations, but these aren't important for my uses as my practical needs are for the SAM UVa lamps that have little output in these regions. ''' for lambda1, weight in weighting_curve: w = Weighting(lambda1) ratio = w/weight print "%d %.2f" % (lambda1, ratio), if 0.8 < ratio < 1.2: print else: print " *" def err(msg): sys.stderr.write(msg) sys.exit(1) def Usage(): print '''Usage: %s op parameters Operates on StellarNet spectrum files. int[egrate] [-w] lambda1 lambda2 file1 [file2...] Prints the integrated irradiance for a StellarNet irradiance file integrated over the indicated wavelengths. If the -w option is given, the effective irradiance over the UV band is printed. dif[ference] file1 file2 output_file Calculates the difference between two spectra and writes it to the output file. ''' % sys.argv[0] sys.exit(1) def ReadLines(file): lines = open(file).readlines() # Remove lines from front of file beginning with '"' and from end of # file beginning with ':'. while lines[0][0] == '"' or lines[0][0] == '#': lines = lines[1:] while lines[-1][0] == ':': lines = lines[0:-1] lines = [s.split() for s in lines] for ix in xrange(len(lines)): lines[ix] = (float(lines[ix][0]), float(lines[ix][1])) assert(len(lines) > 2) return lines def SelectSubset(lines, lambda1, lambda2): '''Return all lines with wavelengths between lambda1 and lambda2 inclusive. ''' return [x for x in lines if lambda1 <= x[0] <= lambda2] def Integrate(Lines, lambda1, lambda2): '''Use the trapezoidal rule. From Bartsch, "Handbook of Mathematical Formulas", Academic Press, 1974, page 361. ''' assert(lambda1 < lambda2) lines = SelectSubset(Lines, lambda1, lambda2) assert(len(lines) > 1) a = lines[0][0] b = lines[-1][0] ya = lines[0][1] yb = lines[-1][1] n = len(lines) - 1 # Number of intervals interior = lines[1:-1] def weight(lambda1): if calculate_weighted: return Weighting(lambda1) else: return 1.0 interior_sum = sum([x[1]*weight(x[0]) for x in interior]) power = (ya + yb + 2*interior_sum) * ((b - a)/(2*n)) return power def PrintResults(lambda1, lambda2, file): power = Integrate(ReadLines(file), lambda1, lambda2) print "%14.6f %s" % (power, file) def Test(): '''We should get 5.7411 for the band [300, 400] nm for the following data: ''' data = [ [298.50, 1.840E-004], [299.00, 5.888E-005], [299.50, 8.308E-005], [300.00, 9.661E-005], [300.50, 8.834E-005], [301.00, 6.585E-005], [301.50, 5.899E-006], [302.00, 0.000E+000], [302.50, 0.000E+000], [303.00, 0.000E+000], [303.50, 2.206E-005], [304.00, 5.637E-005], [304.50, 6.309E-005], [305.00, 2.241E-005], [305.50, 1.310E-004], [306.00, 2.665E-004], [306.50, 3.427E-004], [307.00, 3.719E-004], [307.50, 3.749E-004], [308.00, 3.189E-004], [308.50, 2.610E-004], [309.00, 1.326E-004], [309.50, 7.072E-006], [310.00, 0.000E+000], [310.50, 0.000E+000], [311.00, 0.000E+000], [311.50, 3.437E-005], [312.00, 2.718E-004], [312.50, 5.712E-004], [313.00, 7.039E-004], [313.50, 6.865E-004], [314.00, 5.482E-004], [314.50, 3.919E-004], [315.00, 2.656E-004], [315.50, 1.799E-004], [316.00, 1.537E-004], [316.50, 1.633E-004], [317.00, 1.537E-004], [317.50, 1.630E-004], [318.00, 8.804E-005], [318.50, 4.805E-005], [319.00, 3.060E-005], [319.50, 2.281E-006], [320.00, 1.373E-006], [320.50, 5.559E-005], [321.00, 1.739E-004], [321.50, 3.332E-004], [322.00, 4.217E-004], [322.50, 3.234E-004], [323.00, 3.017E-004], [323.50, 3.125E-004], [324.00, 3.208E-004], [324.50, 4.403E-004], [325.00, 4.846E-004], [325.50, 5.226E-004], [326.00, 5.530E-004], [326.50, 5.528E-004], [327.00, 5.726E-004], [327.50, 6.793E-004], [328.00, 7.918E-004], [328.50, 9.023E-004], [329.00, 9.472E-004], [329.50, 1.030E-003], [330.00, 1.184E-003], [330.50, 1.249E-003], [331.00, 1.325E-003], [331.50, 1.316E-003], [332.00, 1.389E-003], [332.50, 1.599E-003], [333.00, 1.905E-003], [333.50, 2.236E-003], [334.00, 2.581E-003], [334.50, 2.905E-003], [335.00, 3.114E-003], [335.50, 3.244E-003], [336.00, 3.290E-003], [336.50, 3.308E-003], [337.00, 3.351E-003], [337.50, 3.435E-003], [338.00, 3.580E-003], [338.50, 3.784E-003], [339.00, 4.050E-003], [339.50, 4.402E-003], [340.00, 4.688E-003], [340.50, 4.878E-003], [341.00, 5.013E-003], [341.50, 5.173E-003], [342.00, 5.453E-003], [342.50, 5.842E-003], [343.00, 6.195E-003], [343.50, 6.531E-003], [344.00, 6.784E-003], [344.50, 7.021E-003], [345.00, 7.290E-003], [345.50, 7.487E-003], [346.00, 7.778E-003], [346.50, 8.180E-003], [347.00, 8.650E-003], [347.50, 9.202E-003], [348.00, 9.782E-003], [348.50, 1.032E-002], [349.00, 1.093E-002], [349.50, 1.160E-002], [350.00, 1.239E-002], [350.50, 1.331E-002], [351.00, 1.435E-002], [351.50, 1.556E-002], [352.00, 1.699E-002], [352.50, 1.881E-002], [353.00, 2.085E-002], [353.50, 2.313E-002], [354.00, 2.571E-002], [354.50, 2.872E-002], [355.00, 3.236E-002], [355.50, 3.653E-002], [356.00, 4.121E-002], [356.50, 4.643E-002], [357.00, 5.242E-002], [357.50, 5.937E-002], [358.00, 6.705E-002], [358.50, 7.526E-002], [359.00, 8.417E-002], [359.50, 9.391E-002], [360.00, 1.043E-001], [360.50, 1.157E-001], [361.00, 1.276E-001], [361.50, 1.406E-001], [362.00, 1.556E-001], [362.50, 1.716E-001], [363.00, 1.878E-001], [363.50, 2.027E-001], [364.00, 2.160E-001], [364.50, 2.297E-001], [365.00, 2.464E-001], [365.50, 2.644E-001], [366.00, 2.798E-001], [366.50, 2.890E-001], [367.00, 2.919E-001], [367.50, 2.907E-001], [368.00, 2.887E-001], [368.50, 2.868E-001], [369.00, 2.850E-001], [369.50, 2.830E-001], [370.00, 2.808E-001], [370.50, 2.779E-001], [371.00, 2.742E-001], [371.50, 2.698E-001], [372.00, 2.639E-001], [372.50, 2.568E-001], [373.00, 2.487E-001], [373.50, 2.407E-001], [374.00, 2.326E-001], [374.50, 2.243E-001], [375.00, 2.154E-001], [375.50, 2.056E-001], [376.00, 1.954E-001], [376.50, 1.851E-001], [377.00, 1.753E-001], [377.50, 1.662E-001], [378.00, 1.577E-001], [378.50, 1.495E-001], [379.00, 1.413E-001], [379.50, 1.331E-001], [380.00, 1.248E-001], [380.50, 1.168E-001], [381.00, 1.093E-001], [381.50, 1.024E-001], [382.00, 9.622E-002], [382.50, 9.050E-002], [383.00, 8.499E-002], [383.50, 7.968E-002], [384.00, 7.463E-002], [384.50, 6.979E-002], [385.00, 6.524E-002], [385.50, 6.072E-002], [386.00, 5.622E-002], [386.50, 5.195E-002], [387.00, 4.800E-002], [387.50, 4.451E-002], [388.00, 4.151E-002], [388.50, 3.867E-002], [389.00, 3.583E-002], [389.50, 3.303E-002], [390.00, 3.017E-002], [390.50, 2.768E-002], [391.00, 2.548E-002], [391.50, 2.350E-002], [392.00, 2.180E-002], [392.50, 2.016E-002], [393.00, 1.869E-002], [393.50, 1.726E-002], [394.00, 1.584E-002], [394.50, 1.453E-002], [395.00, 1.330E-002], [395.50, 1.210E-002], [396.00, 1.105E-002], [396.50, 1.004E-002], [397.00, 9.152E-003], [397.50, 8.383E-003], [398.00, 7.629E-003], [398.50, 7.000E-003], [399.00, 6.375E-003], [399.50, 5.734E-003], [400.00, 5.214E-003], [400.50, 4.736E-003], [401.00, 4.320E-003], [401.50, 3.992E-003], [402.00, 3.689E-003], ] irradiance = Integrate(data, 300, 400) expected = 5.7411 # W/m^2 assert(abs(irradiance - expected) < 1e-4) print "Test passed" def ProcessCommandLine(): if len(sys.argv) < 2: Usage() command = sys.argv[1][:3] if command == "int": try: optlist, args = getopt.getopt(sys.argv[2:], "w") except getopt.error, str: print "getopt error: %s\n" % str sys.exit(1) for opt in optlist: if opt[0] == "-w": global calculate_weighted calculate_weighted = 1 if len(args) < 3: Usage() try: lambda1 = float(args[0]) lambda2 = float(args[1]) if lambda1 >= lambda2: err("lambda1 must be < lambda2\n") if lambda1 <= 0: err("lambda1 must be > 0\n") if lambda2 <= 0: err("lambda2 must be > 0\n") except: Usage() elif command == "dif": # Need dif file1 file2 output_file if len(sys.argv) < 5: Usage() else: args = sys.argv[2:] else: Usage() args = [command] + args return args def CalculateIrradiance(args): lambda1 = float(args[1]) lambda2 = float(args[2]) print "Irradiance in W/m^2 over [%.1f, %.1f] nm for the following files:" % \ (lambda1, lambda2) for file in args[3:]: PrintResults(lambda1, lambda2, file) def CalculateDifference(args): # Put spectra into dictionary keyed by wavelength wavelength1 = {} wavelength2 = {} for wavelength, signal in ReadLines(args[1]): wavelength1[wavelength] = signal for wavelength, signal in ReadLines(args[2]): wavelength2[wavelength] = signal # Calculate difference diff = {} for wavelength in wavelength1: if wavelength2.has_key(wavelength): diff[wavelength] = wavelength1[wavelength] - wavelength2[wavelength] else: diff[wavelength] = wavelength1[wavelength] # Change difference dictionary back into a list list = diff.keys() list.sort() for ix in xrange(len(list)): list[ix] = (list[ix], diff[list[ix]]) # Write to output file ofp = open(args[3], "wb") for item in list: ofp.write("%.1f %.4e\n" % item) def main(): args = ProcessCommandLine() if args[0] == "int": CalculateIrradiance(args) elif args[0] == "dif": CalculateDifference(args) main()