' --------------------------------------------------------------------------------------------
' polarizationAnalysis.frs
' This file contains all of the FRED routines used to perform polarization analysis for VIIRS.
' Many of the subroutines in this file were written to support specific analysis requests and
' are not commonly used. The most commonly used routines are as follows (given in the order in
' which they are most commonly called; first routine calls the second, which calls the third,
' etc.):
' 	Main: defines the source spectrum, waveband to analyze, and scan angles for the analysis.
'      Generally, this is the highest level routine used, and as such requires user input.
'	analyzeAtSevenScanMirrorAngles (or similar): loops over the scan mirror angles
'   computePolarizationStatisticsForBand (or similar): sets the waveband, calls
'      setupCoatingsAndSources, loops over each detector in the band,
'      calls computePolarizationStatisticsSixSources for each one
'   setupCoatingsAndSources: a very important routine. assigns coatings to each surface,
'      which means that it defines the system (F1, J1, J2) being analyzed. All coating
'      assignments happen in this routine. Also sets the appropriate sources for the
'      FPA being analyzed (VISNIR, SWMWIR, LWIR).
'   computePolarizationStatisticsSixSources (or similar): does the ray trace, computes the
'      polarization sensitivity.
'
' Questions about the model and script can be sent to Eric Fest, eric@raytheon.com.
' Please note any major changes to this file in the history below.
' --------------------------------------------------------------------------------------------
' Copyright (c) Raytheon Proprietary 2014
' This FRED model and script may not be used without permission from Raytheon Company.
' --------------------------------------------------------------------------------------------
' 20090420: Initial released version (ECF)
' 20130125: Ported to WinWrap BASIC (ECF)
' 20140430: Added the above header info (ECF)
' 20140801: Modified Eric Fest's code to only look at polarized transmittance for combinations desired.
' --------------------------------------------------------------------------------------------

Option Explicit
' The const variables below are used in the call to setSourceWvlAndDirection and determine which sources get set to traceable.
Const unpolarizedFullPupil = 0			' used for system transmittance studies
Const sixPolarizationsFullPupil = 1		' used for DoLP studies. This is the minimum number of sources required to compute the Stokes vector.
Const onePolarizationSubPupil = 2		' used for DoLP studies to compute the pupil dependency of the polarization.

Enum spectrumType
	List
	Spectra
End Enum

Type spectrumData	' every piece of data needed to set the source spectrum
	wvbandStr As String	' M1-M7, and I1
	spectTyp As spectrumType	' Either 'List' or 'Spectra'
	wavelengths() As Double	' the wavelengths for "List" type. Also, its size is number of wavelengths for a "Spectra" Type
	spectID As Long	' spectrum ID for a "Spectra" type
End Type

Sub main
	' These are all the global variables we need to perform any analysis
	Dim covar(3,3) As Double, alpha(3,3) As Double, oneda(3,1) As Double	' For curve fitting
	Dim skipDOCP As Boolean
	Dim D2DDOLP(16) As Double, D2DModulation(48) As Double, DOLP As Double, modulation(4) As Double
	Dim maxFrpVariation As Double, maxFPpVariation As Double
	Dim xi(1 To 54, 1 To 54) As Double	' for powell's method
	Dim ncom As Integer ' for powell's method
	Dim pcom() As Double, xicom() As Double ' for powell's method

	' Global variables for simplex optimization
	Dim VertexMat(55,54) As Double       ' Amoeba Vertices (NDim+1, NDim)
    	                               ' Size has to be (NDim+1, NDim)
	Dim IdentityMatrix(55,54) As Integer ' Identity matrix
                                   ' Has to be same size as VertexMax
	Dim Variables() As Double          ' Array to store initial design, used in optimization
	Dim Increments() As Double         ' Array to store initial variable increments
	Dim availableVariables() As Double ' All available variables
	Dim maxAvailableVariables() As Double	' Max value of available variables
	Dim activeAvailableVariables() As Boolean   ' Which available variables are being used by this merit function
	Dim NDim As Integer                ' Number of dimensions (variables)
	Dim MaxLoop As Integer             ' Most number of iterations
	Dim VolumeTolerance As Double      ' Largest space enclosed by vertices
	Dim optimizationOutputFilename As String, muellerMatrixFilename As String, detBaseStr As String, detFullStr As String
	Dim muellerMatrixFilenameRoot As String
	Dim pathID As Long, detNo As Long, detID As Long, nRaysAcrossPupil As Long
	Dim advr As T_ADVANCEDRAYTRACE
	Dim collectorID As Long, i As Long
	Dim detEllipticity As Double, detOrientation As Double
	Dim detRightHanded As Boolean
	Dim advRaytrace As T_ADVANCEDRAYTRACE
	Dim thinFilmLayerThicknessScaleFactorsForOpt(10, 2)	As Double ' Thickness scale factors for thin film optimization for (surface #, layer #). New thickness = Old thickness * (1 +/- scale factor)
	Dim maxAbsThicknessScaleFactor As Double
	Dim anType As Long
	Dim spectData As spectrumData
	Dim lowBoundWv As Double, highBoundWv As Double, delWv As Double

	'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	spectData.spectTyp = Spectra
	ReDim spectData.wavelengths(25)
	' Specify spectrum to use
	spectData.spectID = SpectrumFindName( "Flat" )
	'Make sure you're using the HAM coating you want
	'Right now the HAM A vs HAM B difference is turning out to be a wash of +/-1% error
	'A good portion of the data is in the excel if you feel like digging through a million plots
	'wvbandStr = "M4"
	'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	' RSR
	'spectData.spectTyp = List
	'ReDim spectData.wavelengths(25)

	'spectData.wvbandStr = "M1"
	'lowBoundWv = 0.3928 '- edgeBuffer
	'highBoundWv = 0.4236 '+ edgeBuffer
	'spectData.wvbandStr = "M2"
	'lowBoundWv = 0.43 '- edgeBuffer
	'highBoundWv = 0.4582 '+ edgeBuffer
	'spectData.wvbandStr = "M3"
	'lowBoundWv = 0.4718 '- edgeBuffer
	'highBoundWv = 0.5032 '+ edgeBuffer
	'spectData.wvbandStr = "M5"
	'lowBoundWv = 0.5384 '- edgeBuffer
	'highBoundWv = 0.571 '+ edgeBuffer
	'spectData.wvbandStr = "M5"
	'lowBoundWv = 0.654 '- edgeBuffer
	'highBoundWv = 0.687 '+ edgeBuffer
	'spectData.wvbandStr = "M6"
	'lowBoundWv = 0.7322 '- edgeBuffer
	'highBoundWv = 0.7574 '+ edgeBuffer
	'spectData.wvbandStr = "M7"
	'lowBoundWv = 0.8376 '- edgeBuffer
	'highBoundWv = 0.8892 '+ edgeBuffer
	'spectData.wvbandStr = "I1"
	'lowBoundWv =  0.5914 '- edgeBuffer
	'highBoundWv = 0.6874 '+ edgeBuffer
	'spectData.wvbandStr = "I2"
	'lowBoundWv = 0.8376 '- edgeBuffer
	'highBoundWv = 0.8892 '+ edgeBuffer

	delWv = (highBoundWv - lowBoundWv) / 25
	For i = 0 To 24
		spectData.wavelengths(i) = lowBoundWv + delWv*0.5 + i*delWv
	Next i
	Call computeRelativeSpectralResponse(spectData)

	'Call computePolarizationEvolutionSixSources()
End Sub

Sub computeRelativeSpectralResponse(spectData As spectrumData)
	Print "Scan angle = -55"
	Call setRTAScanAngle(-101.)	' -55
	Call computeRelativeSpectralResponseForBand(spectData)
	Print "Scan angle = -45"
	Call setRTAScanAngle(-91.)	' -45
	Call computeRelativeSpectralResponseForBand(spectData)
	Print "Scan angle = -37"
	Call setRTAScanAngle(-83.)	' -37
	Call computeRelativeSpectralResponseForBand(spectData)
	Print "Scan angle = -30"
	Call setRTAScanAngle(-76.)	' -30
	Call computeRelativeSpectralResponseForBand(spectData)
	Print "Scan angle = -20"
	Call setRTAScanAngle(-66.)	' -20
	Call computeRelativeSpectralResponseForBand(spectData)
	Print "Scan angle = -15"
	Call setRTAScanAngle(-61.)	' -15
	Call computeRelativeSpectralResponseForBand(spectData)
	Print "Scan angle = -8"
	Call setRTAScanAngle(-54.)	' -8
	Call computeRelativeSpectralResponseForBand(spectData)
	Print "Scan angle = 4"
	Call setRTAScanAngle(-42.)	' +4
	Call computeRelativeSpectralResponseForBand(spectData)
	Print "Scan angle = 22"
	Call setRTAScanAngle(-24.)	' +22
	Call computeRelativeSpectralResponseForBand(spectData)
	Print "Scan angle = 45"
	Call setRTAScanAngle(-1.)	' +45
	Call computeRelativeSpectralResponseForBand(spectData)
	Print "Scan angle = 55"
	Call setRTAScanAngle(9.)	' +55
	Call computeRelativeSpectralResponseForBand(spectData)
End Sub

Sub computeRelativeSpectralResponseForBand(spectData As spectrumData)
	Dim detNo() As Long, nRays As Long, pathID As Long, i As Long, j As Long, detID As Long
	Dim skipDOCP As Boolean
	Dim spectDataOneWavelength As spectrumData
	Dim DOLP As Double, modulation(4) As Double
	Dim detBaseStr As String, detFullStr As String

	' Define the detectors of interest
	If spectData.wvbandStr = "I1" Or spectData.wvbandStr = "I2" Then
	   ReDim detNo(16)
	   detNo(0) = 1
	   detNo(1) = 2
	   detNo(2) = 3
	   detNo(3) = 4
	   detNo(4) = 5
	   detNo(5) = 6
	   detNo(6) = 7
	   detNo(7) = 8
	   detNo(8) = 9
	   detNo(9) = 10
	   detNo(10) = 11
	   detNo(11) = 12
	   detNo(12) = 13
	   detNo(13) = 14
	   detNo(14) = 15
	   detNo(15) = 16
	   detNo(16) = 17
	   detNo(17) = 18
	   detNo(18) = 19
	   detNo(19) = 20
	   detNo(20) = 21
	   detNo(21) = 22
	   detNo(22) = 23
	   detNo(23) = 24
	   detNo(24) = 25
	   detNo(25) = 26
	   detNo(26) = 27
	   detNo(27) = 28
	   detNo(28) = 29
	   detNo(29) = 30
	   detNo(30) = 31
	   detNo(31) = 32
	Else
		ReDim detNo(16)
		detNo(0) = 1
		detNo(1) = 2
		detNo(2) = 3
		detNo(3) = 4
		detNo(4) = 5
		detNo(5) = 6
		detNo(6) = 7
		detNo(7) = 8
		detNo(8) = 9
		detNo(9) = 10
		detNo(10) = 11
		detNo(11) = 12
		detNo(12) = 13
		detNo(13) = 14
		detNo(14) = 15
		detNo(15) = 16
	End If

	' Initialize
	nRays = 21
	spectDataOneWavelength.wvbandStr = spectData.wvbandStr
	spectDataOneWavelength.spectTyp = List
	ReDim spectDataOneWavelength.wavelengths(1)
	detBaseStr = "Geometry.VISNIR_detectors." + spectData.wvbandStr
	detBaseStr = detBaseStr + "_detectors.det"
	Print "Band = "; spectData.wvbandStr
	Print "Wavelength (um)" & Chr(9) & "det num" & Chr(9) & "P0" & Chr(9) & "P45" & Chr(9) & "P90" & Chr(9) & "P135"

	' Loop over the detectors
	For i = 0 To UBound(detNo) - 1

		detFullStr = detBaseStr + CStr(detNo(i))
		detID = FindFullName(detFullStr)

		' Loop over the wavelengths
		For j = 0 To UBound(spectData.wavelengths) - 1

			spectDataOneWavelength.wavelengths(0) = spectData.wavelengths(j)
			'Print "Wavelength = "; spectDataOneWavelength.wavelengths(0)

			' Setup the coatings and sources for this wavelength and det num
			Call setupCoatingsAndSources(spectDataOneWavelength, nRays, detNo(i)-1, skipDOCP, pathID)

			' Do the ray trace
			Call computePolarizationStatisticsTwelveSources(detID, pathID, DOLP, modulation)

		Next j
	Next i
End Sub

' This script sets the wavelength and direction cosines of the Collimated Test Beam sources for the VIIRS polarization model.
' The direction is defined in terms of alpha (rotation about x-axis, in track dir) and beta (rotation about y-axis, in scan dir).
'Sub setupSixSources(waveband As String, numRays As Long, alpha1 As Double, skipDOCP As Boolean, pathID As Long)
Sub setupCoatingsAndSources(spectData As spectrumData, numRays As Long, detNo As Long, skipDOCP As Boolean, pathID As Long)

	Dim beta As Double, sourceDiameter As Double, sourceBetaMBands(16) As Double, alpha As Double, sourceBetaIBands(32) As Double
	Dim source0ID As Long, source45ID As Long, source90ID As Long, source135ID As Long, sourceRID As Long, sourceLID As Long, mirrorCoatingID As Long
	Dim DBS1CoatingID As Long, DBS1RaytraceCtrlID As Long, DBS2RaytraceCtrlID As Long, degradedMirrorCoatingID As Long
	Dim RTA1CoatingID As Long, RTA2CoatingID As Long, RTA3CoatingID As Long, FM1CoatingID As Long, HAMCoatingID As Long, FM2CoatingID As Long
	Dim AFT1CoatingID As Long, AFT2CoatingID As Long, AFT3CoatingID As Long, AFT4CoatingID As Long, i As Long
	Dim source15ID As Long, source30ID As Long, source60ID As Long, source75ID As Long, source105ID As Long, source120ID As Long, source150ID As Long, source165ID As Long
	Dim numWv As Long, numWvForSource As Long, spectWeight As String
	Dim skipSetCoatings As Boolean
	Dim sourceMMID As Long, surfID As Long, onWvPerBandVALUE As String
	Dim edgeBuffer As Double, lowBoundWv As Double, highBoundWv As Double
	Dim M1CoatingID As Long, M2CoatingID As Long, M3CoatingID As Long, M4CoatingID As Long, M5CoatingID As Long
	Dim M6CoatingID As Long, M7CoatingID As Long, I1CoatingID As Long, I2CoatingID As Long
	Dim angleTuneUp As Double
	Dim wvSpec As T_WAVELENGTHSPEC

	skipSetCoatings = False
	angleTuneUp = 0.
	If spectData.wvbandStr = "I1" Then
		sourceBetaIBands(0)  = -.39783  : sourceBetaIBands(1)  = -.37216  : sourceBetaIBands(2)  =-.3465  : sourceBetaIBands(3)  =-.32083
		sourceBetaIBands(4)  = -.29517  : sourceBetaIBands(5) = -.2695    : sourceBetaIBands(6) =-.24383  : sourceBetaIBands(7)  =-.21817
		sourceBetaIBands(8)  = -.1925   : sourceBetaIBands(9)  = -.16683  : sourceBetaIBands(10) =-.14116 : sourceBetaIBands(11) =-.1155
		sourceBetaIBands(12) = -.08983  : sourceBetaIBands(13) = -.06416  : sourceBetaIBands(14) =-.0385  : sourceBetaIBands(15) =-.01283
		sourceBetaIBands(16)  = .01283  : sourceBetaIBands(17)  = .0385   : sourceBetaIBands(18)  =.06416 : sourceBetaIBands(19) = .08983
		sourceBetaIBands(20)  = .1155   : sourceBetaIBands(21) = .14116  : sourceBetaIBands(22) = .16683  : sourceBetaIBands(23) = .1925
		sourceBetaIBands(24)  = .21817  : sourceBetaIBands(25)  = .24383  : sourceBetaIBands(26) =.2695   : sourceBetaIBands(27) = .29517
		sourceBetaIBands(28) = .32083   : sourceBetaIBands(29) = .3465    : sourceBetaIBands(30) =.37216  : sourceBetaIBands(31) = .39783
		beta = sourceBetaIBands(detNo)
	Else
		sourceBetaMBands(0)  = -.385    : sourceBetaMBands(1)  = -.33367 : sourceBetaMBands(2)  =-.28233 : sourceBetaMBands(3)  =-.231
		sourceBetaMBands(4)  = -.17967  : sourceBetaMBands(5)  = -.12834 : sourceBetaMBands(6)  =-.077   : sourceBetaMBands(7)  =-.025664
		sourceBetaMBands(8)  = .025664  : sourceBetaMBands(9)  = .077    : sourceBetaMBands(10) =.12834  : sourceBetaMBands(11) = .17967
		sourceBetaMBands(12) = .231     : sourceBetaMBands(13) = .28233  : sourceBetaMBands(14) =.33367  : sourceBetaMBands(15) = .385
		beta = sourceBetaMBands(detNo)
	End If

	' First, set up the coatings and raytrace controls so we trace to the appropriate FPA assembly
	spectWeight = spectData.wvbandStr
	If spectData.wvbandStr = "M1" Or spectData.wvbandStr = "M2" Or spectData.wvbandStr = "M3" Or spectData.wvbandStr = "M4" Or _
		spectData.wvbandStr = "M5" Or spectData.wvbandStr = "M6" Or spectData.wvbandStr= "M7" Or _
		spectData.wvbandStr = "I1" Or spectData.wvbandStr = "I2" Then
		' VISNIR
		' Get the coating IDs of the mirrors in the visible
		'RTA1CoatingID = FindCoating("MODIS FSS99-500")	' Thin-film model
		'RTA1CoatingID = FindCoating("MODIS FSS99-500 EOL")	' Thin-film model EOL
		'RTA1CoatingID = FindCoating("F1 RTA VIS (AOI Extrapolated)") ' Nominal
		RTA1CoatingID = FindCoating("J1 FMA VIS (AOI Extrapolated)")	' Nominal
		'RTA1CoatingID = FindCoating("F1 RTA1 VIS Degraded")
		'RTA1CoatingID = FindCoating("MODIS FSS99-500 RTA1")	' For coating optimization
		'RTA1CoatingID = FindCoating("Specification for Mirrors")	' For spec study

		'RTA2CoatingID = FindCoating("MODIS FSS99-500") ' Thin-film model
		'RTA2CoatingID = FindCoating("MODIS FSS99-500 EOL")	' Thin-film model EOL
		'RTA2CoatingID = FindCoating("F1 RTA VIS (AOI Extrapolated)")	' Nominal
		RTA2CoatingID = FindCoating("J1 FMA VIS (AOI Extrapolated)")	' Nominal
		'RTA2CoatingID = FindCoating("MODIS FSS99-500 RTA2")	' For coating optimization
		'RTA2CoatingID = FindCoating("Specification for Mirrors")	' For spec study

		'RTA3CoatingID = FindCoating("MODIS FSS99-500") ' Thin-film model
		'RTA3CoatingID = FindCoating("MODIS FSS99-500 EOL")	' Thin-film model EOL
		'RTA3CoatingID = FindCoating("F1 RTA VIS (AOI Extrapolated)")	' Nominal
		RTA3CoatingID = FindCoating("J1 FMA VIS (AOI Extrapolated)")	' Nominal
		'RTA3CoatingID = FindCoating("F1 RTA3 VIS Degraded")
		'RTA3CoatingID = FindCoating("MODIS FSS99-500 RTA3")	' For coating optimization
		'RTA3CoatingID = FindCoating("Specification for Mirrors")	' For spec study

		'FM1CoatingID = FindCoating("MODIS FSS99-500") ' Thin-film model
		'FM1CoatingID = FindCoating("MODIS FSS99-500 EOL")	' Thin-film model EOL
		'FM1CoatingID = FindCoating("F1 RTA VIS (AOI Extrapolated)")	' Nominal
		FM1CoatingID = FindCoating("J1 FMA VIS (AOI Extrapolated)")	' Nominal
		'FM1CoatingID = FindCoating("MODIS FSS99-500 FM1")	' For coating optimization
		'FM1CoatingID = FindCoating("Specification for Mirrors")	' For spec study

		'HAMCoatingID = FindCoating("MODIS FSS99-500") ' Thin-film model
		'HAMCoatingID = FindCoating("MODIS FSS99-500 EOL (20 nm Kapton)")	' Thin-film model EOL 1
		'HAMCoatingID = FindCoating("MODIS FSS99-500 EOL (20 nm Polyethylene)")	' Thin-film model EOL 2
		'HAMCoatingID = FindCoating("MODIS FSS99-500 EOL (5 nm Polyethylene)")	' Thin-film model EOL 3
		'HAMCoatingID = FindCoating("MODIS FSS99-500 EOL (10 nm Polyethylene)")	' Thin-film model EOL 4
		'HAMCoatingID = FindCoating("F1 HAM VIS Corrected")	' Nominal
		HAMCoatingID = FindCoating("J1 FMA VIS (AOI Extrapolated)")	' Nominal
		'HAMCoatingID = FindCoating( "J1 HAM A VIS (AOI Extrapolated)" )
		'HAMCoatingID = FindCoating( "J1 HAM B VIS (AOI Extrapolated)" )
		'HAMCoatingID = FindCoating("F1 HAM VIS Degraded")
		'HAMCoatingID = FindCoating("MODIS FSS99-500 HAM")	' For coating optimization
		'HAMCoatingID = FindCoating("Specification for Mirrors")	' For spec study

		'FM2CoatingID = FindCoating("MODIS FSS99-500") ' Thin-film model
		'FM2CoatingID = FindCoating("MODIS FSS99-500 EOL")	' Thin-film model EOL
		'FM2CoatingID = FindCoating("F1 RTA VIS (AOI Extrapolated)")	' Nominal
		FM2CoatingID = FindCoating("J1 FMA VIS (AOI Extrapolated)")	' Nominal
		'FM2CoatingID = FindCoating("F1 FM2 VIS Degraded")
		'FM2CoatingID = FindCoating("MODIS FSS99-500 FM2")	' For coating optimization
		'FM2CoatingID = FindCoating("Specification for Mirrors")	' For spec study

		'AFT1CoatingID = FindCoating("MODIS FSS99-500") ' Thin-film model
		'AFT1CoatingID = FindCoating("MODIS FSS99-500 EOL")	' Thin-film model EOL
		'AFT1CoatingID = FindCoating("F1 FMA VIS (AOI Extrapolated)")	' Nominal
		AFT1CoatingID = FindCoating("J1 FMA VIS (AOI Extrapolated)")	' Nominal
		'AFT1CoatingID = FindCoating("F1 AFT1 VIS Degraded")
		'AFT1CoatingID = FindCoating("MODIS FSS99-500 AFT1")	' For coating optimization
		'AFT1CoatingID = FindCoating("Specification for Mirrors")	' For spec study

		'AFT2CoatingID = FindCoating("MODIS FSS99-500") ' Thin-film model
		'AFT2CoatingID = FindCoating("MODIS FSS99-500 EOL")	' Thin-film model EOL
		'AFT2CoatingID = FindCoating("F1 FMA VIS (AOI Extrapolated)")	' Nominal
		AFT2CoatingID = FindCoating("J1 FMA VIS (AOI Extrapolated)")	' Nominal
		'AFT2CoatingID = FindCoating("F1 AFT2 VIS Degraded")
		'AFT2CoatingID = FindCoating("Absorb")
		'AFT2CoatingID = FindCoating("MODIS FSS99-500 AFT2")	' For coating optimization
		'AFT2CoatingID = FindCoating("Specification for Mirrors")	' For spec study

		'AFT3CoatingID = FindCoating("MODIS FSS99-500") ' Thin-film model
		'AFT3CoatingID = FindCoating("MODIS FSS99-500 EOL")	' Thin-film model EOL
		'AFT3CoatingID = FindCoating("F1 FMA VIS (AOI Extrapolated)")	' Nominal
		AFT3CoatingID = FindCoating("J1 FMA VIS (AOI Extrapolated)")	' Nominal
		'AFT3CoatingID = FindCoating("F1 AFT3 VIS Degraded")
		'AFT3CoatingID = FindCoating("MODIS FSS99-500 AFT3")	' For coating optimization
		'AFT3CoatingID = FindCoating("Specification for Mirrors")	' For spec study

		'AFT4CoatingID = FindCoating("MODIS FSS99-500") ' Thin-film model
		'AFT4CoatingID = FindCoating("MODIS FSS99-500 EOL")	' Thin-film model EOL
		'AFT4CoatingID = FindCoating("F1 FMA VIS (AOI Extrapolated)")	' Nominal
		AFT4CoatingID = FindCoating("J1 FMA VIS (AOI Extrapolated)")	' Nominal
		'AFT4CoatingID = FindCoating("F1 AFT4 VIS Degraded")
		'AFT4CoatingID = FindCoating("Absorb")
		'AFT4CoatingID = FindCoating("MODIS FSS99-500 AFT4")	' For coating optimization
		'AFT4CoatingID = FindCoating("Specification for Mirrors")	' For spec study

		' Get the coating ID of the DBS1 data in the visible
		'DBS1CoatingID = FindCoating("F1 DBS1 VIS Predicted")	' Nominal
		'DBS1CoatingID = FindCoating("F1 DBS1 VIS Corrected")
		DBS1CoatingID = FindCoating("J1 DBS1 VIS Predicted")	' Nominal
		'DBS1CoatingID = FindCoating( "J1 DBS1 with half (rs-rp)" )
		'DBS1CoatingID = FindCoating("F1 DBS1 VIS Degraded")
		'DBS1CoatingID = FindCoating("Specification for DBS1")	' For spec study

		' Get the reflect specular raytrace ctrl ID for DBS1
		DBS1RaytraceCtrlID = FindRaytraceCtrl("Reflect Specular")

		' Set the IFA coatings

		'M1CoatingID = FindCoating("F1 M1")
		'M2CoatingID = FindCoating("F1 M2")
		'M3CoatingID = FindCoating("F1 M3")
		'M4CoatingID = FindCoating("F1 M4")
		'M5CoatingID = FindCoating("F1 M5")
		'M6CoatingID = FindCoating("F1 M6")
		'M7CoatingID = FindCoating("F1 M7")
		'I1CoatingID = FindCoating("F1 I1")

		M1CoatingID = FindCoating("J1 M1")
		M2CoatingID = FindCoating("J1 M2")
		M3CoatingID = FindCoating("J1 M3")
		M4CoatingID = FindCoating("J1 M4")
		M5CoatingID = FindCoating("J1 M5")
		M6CoatingID = FindCoating("J1 M6")
		M7CoatingID = FindCoating("J1 M7-I2")
		I1CoatingID = FindCoating("J1 I1")

		'M1CoatingID = FindCoating("J2 M1")
		'M1CoatingID = FindCoating("J2 M1 New") 'Materion New J2 M1 coating
		'M1CoatingID = FindCoating( "J2 M1 New X 0.5" ) ' Vijays half-polarization coating
		'M2CoatingID = FindCoating("J2 M2")
		'M2CoatingID = FindCoating("J2 M2 New") 'Materion New J2 M2 Coating
		'M3CoatingID = FindCoating("J2 M3")
		'M3CoatingID = FindCoating("J2 M3 New") 'Materion New J2 M3 Coating 9Apr2014
		'M4CoatingID = FindCoating("J2 M4")
		'M4CoatingID = FindCoating("J2 M4 New") 'Materion New J2 M4 Coating 9Apr2014
		'M5CoatingID = FindCoating("J2 M5")
		'M6CoatingID = FindCoating("J2 M6")
		'M7CoatingID = FindCoating("J2 M7-I2")
		'I1CoatingID = FindCoating("J2 I1")

		' Doesn't matter what DBS2 does
		DBS2RaytraceCtrlID = FindRaytraceCtrl("Reflect Specular")

	ElseIf spectData.wvbandStr = "M8" Or spectData.wvbandStr = "M9" Or spectData.wvbandStr = "M10" Or _
		spectData.wvbandStr= "M11" Or spectData.wvbandStr = "M12" Or spectData.wvbandStr = "M13" Or _
		spectData.wvbandStr = "I3" Or spectData.wvbandStr = "I4" Then
		' SWMWIR
		' Get the coating ID of the mirror data in the IR
		mirrorCoatingID = FindCoating("FSS99-500 IR")

		' Get the coating ID of the DBS1 data in the IR
		DBS1CoatingID = FindCoating("Dichroic1 entrance")	' This is all we've got for now

		' Get the transmit specular raytrace ctrl ID for DBS1
		DBS1RaytraceCtrlID = FindRaytraceCtrl("Transmit Specular")

		' Get the reflect specular raytrace ctrl ID for DBS1
		DBS2RaytraceCtrlID = FindRaytraceCtrl("Reflect Specular")

	ElseIf spectData.wvbandStr = "M14" Or spectData.wvbandStr = "M15" Or spectData.wvbandStr = "M16" Or spectData.wvbandStr = "I5" Then
		' Get the coating ID of the mirror data in the IR
		mirrorCoatingID = FindCoating("FSS99-500 IR")

		' Get the coating ID of the DBS1 data in the IR
		DBS1CoatingID = FindCoating("Dichroic1 entrance")	' This is all we've got for now

		' Get the transmit specular raytrace ctrl ID for DBS1
		DBS1RaytraceCtrlID = FindRaytraceCtrl("Transmit Specular")

		' Get the transmit specular raytrace ctrl ID for DBS1
		DBS2RaytraceCtrlID = FindRaytraceCtrl("Transmit Specular")

	End If

	' Set the coatings
	If Not skipSetCoatings Then
		' Primary
		surfID = FindFullName("Geometry.RTA.Fore M1.Surf 6")
		SetSurfCoating surfID, RTA1CoatingID
		' Secondary
		surfID = FindFullName("Geometry.RTA.Fore M2.Surf 8")
		SetSurfCoating surfID, RTA2CoatingID
		' Tertiary
		surfID = FindFullName("Geometry.RTA.Fore M3.Surf 11")
		SetSurfCoating surfID, RTA3CoatingID
		' Fold 1
		surfID = FindFullName("Geometry.RTA.Fold 1.Surf 14")
		SetSurfCoating surfID, FM1CoatingID
		' HAM
		surfID = FindFullName("Geometry.HAM.Ham.Surf 17")
		SetSurfCoating surfID, HAMCoatingID
		' Fold 2
		surfID = FindFullName("Geometry.Fixed Relay Optics.Fold 2.Surf 22")
		SetSurfCoating surfID, FM2CoatingID
		'SetSurfCoating surfID, degradedMirrorCoatingID
		' Aft M1
		surfID = FindFullName("Geometry.Fixed Relay Optics.Aft M1.Surf 25")
		SetSurfCoating surfID, AFT1CoatingID
		'SetSurfCoating surfID, degradedMirrorCoatingID
		' Aft M2
		surfID = FindFullName("Geometry.Fixed Relay Optics.Aft M2.Surf 27")
		SetSurfCoating surfID, AFT2CoatingID
		' Aft M3
		surfID = FindFullName("Geometry.Fixed Relay Optics.Aft M3.Surf 30")
		SetSurfCoating surfID, AFT3CoatingID
		' Aft M4
		surfID = FindFullName("Geometry.Fixed Relay Optics.Aft M4.Surf 32")
		SetSurfCoating surfID, AFT4CoatingID
		' DBS1
		surfID = FindFullName("Geometry.swir.D1.Surf 36")
		SetSurfCoating surfID, DBS1CoatingID
		' IFA
		surfID = FindFullName( "Geometry.VISNIR_Integrated_Filter_Assembly.M1.Entrance Face" )
		SetSurfCoating surfID, M1CoatingID
		surfID = FindFullName( "Geometry.VISNIR_Integrated_Filter_Assembly.M2.Entrance Face" )
		SetSurfCoating surfID, M2CoatingID
		surfID = FindFullName( "Geometry.VISNIR_Integrated_Filter_Assembly.M3.Entrance Face" )
		SetSurfCoating surfID, M3CoatingID
		surfID = FindFullName( "Geometry.VISNIR_Integrated_Filter_Assembly.M4.Entrance Face" )
		SetSurfCoating surfID, M4CoatingID
		surfID = FindFullName( "Geometry.VISNIR_Integrated_Filter_Assembly.M5.Entrance Face" )
		SetSurfCoating surfID, M5CoatingID
		surfID = FindFullName( "Geometry.VISNIR_Integrated_Filter_Assembly.M6.Entrance Face" )
		SetSurfCoating surfID, M6CoatingID
		surfID = FindFullName( "Geometry.VISNIR_Integrated_Filter_Assembly.I2/M7.Entrance Face" )
		SetSurfCoating surfID, M7CoatingID
		surfID = FindFullName( "Geometry.VISNIR_Integrated_Filter_Assembly.I1.Entrance Face" )
		SetSurfCoating surfID, I1CoatingID

		' Set the raytrace controls
		' DBS1
		surfID = FindFullName("Geometry.swir.D1.Surf 36")
		SetSurfRaytraceCtrl surfID, DBS1RaytraceCtrlID
		' DBS2
		surfID = FindFullName("Geometry.lwir.D2.Surf 41")
		SetSurfRaytraceCtrl surfID, DBS2RaytraceCtrlID
	End If

	' Setup the sources
	SetUserDefinedRayPathMode True

	' Set the source angle, Path ID, and ray color
	If spectData.wvbandStr = "M1" Then
		alpha = .608	' About at the center of M1, det 8 (w/alpha=0)
		pathID = PathFind("VISNIR M1")
		wvSpec.colorR = 0
		wvSpec.colorG = 0
		wvSpec.colorB = 255
	ElseIf spectData.wvbandStr = "M2" Then
		alpha = .555	' About at the center of M2, det 8 (w/alpha=0)
		pathID = PathFind("VISNIR M2")
		wvSpec.colorR = 0
		wvSpec.colorG = 0
		wvSpec.colorB = 255
	ElseIf spectData.wvbandStr = "M3" Then
		alpha = .449	' About at the center of M3, det 8 (w/alpha=0)
		pathID = PathFind("VISNIR M3")
		wvSpec.colorR = 0
		wvSpec.colorG = 255
		wvSpec.colorB = 0
	ElseIf spectData.wvbandStr = "M4" Then
		alpha = .5	' About at the center of M4, det 8 (w/alpha=0)
		pathID = PathFind("VISNIR M4")
		wvSpec.colorR = 0
		wvSpec.colorG = 255
		wvSpec.colorB = 0
	ElseIf spectData.wvbandStr = "M5" Then
		alpha = .29	' About at the center of M5, det 8 (w/alpha=0)
		pathID = PathFind("VISNIR M5")
		wvSpec.colorR = 0
		wvSpec.colorG = 255
		wvSpec.colorB = 0
	ElseIf spectData.wvbandStr = "M6" Then
		alpha = .23	' About at the center of M6, det 8 (w/alpha=0)
		pathID = PathFind("VISNIR M6")
		wvSpec.colorR = 255
		wvSpec.colorG = 0
		wvSpec.colorB = 0
	ElseIf spectData.wvbandStr = "M7" Then
		alpha = .34
		pathID = PathFind("VISNIR M7")
		wvSpec.colorR = 255
		wvSpec.colorG = 0
		wvSpec.colorB = 0
	ElseIf spectData.wvbandStr = "I1" Then
		alpha = .4	' About at the center of I1, det 8 (w/alpha=0)
		pathID = PathFind("VISNIR I1")
		wvSpec.colorR = 0
		wvSpec.colorG = 0
		wvSpec.colorB = 255
	ElseIf spectData.wvbandStr = "I2" Then
		'beta = 0
		pathID = PathFind("VISNIR I2")
		wvSpec.colorR = 0
		wvSpec.colorG = 0
		wvSpec.colorB = 255
	ElseIf spectData.wvbandStr = "I1" Then
		alpha = .4	' About at the center of I1, det 8 (w/alpha=0)
		wvSpec.colorR = 255
		wvSpec.colorG = 0
		wvSpec.colorB = 0
	ElseIf spectData.wvbandStr = "I2" Then
		'beta = 0
		wvSpec.colorR = 255
		wvSpec.colorG = 0
		wvSpec.colorB = 0
	End If

	If spectData.spectTyp = Spectra Then

		' We add an "edge buffer" to test for repeatability with the other method of specifying a list of weighted wavelengths.
		' There's a slight, documented difference how each method samples the range of wavelengths, and this should make
		' spectra results do the same sampling as evenly spaced weights with a list.
		If spectData.wvbandStr = "M1" Then
			'edgeBuffer = (0.4236 - 0.3928)/(numWvPerBand*2)
			lowboundWv = 0.3928 '- edgeBuffer
			highboundWv = 0.4236 '+ edgeBuffer
		ElseIf spectData.wvbandStr = "M2" Then
			'edgeBuffer = (0.4582 - 0.43)/(numWvPerBand*2)
			lowboundWv = 0.43 '- edgeBuffer
			highboundWv = 0.4582 '+ edgeBuffer
		ElseIf spectData.wvbandStr = "M3" Then
			'edgeBuffer = (0.5032 - 0.4718)/(numWvPerBand*2)
			lowboundWv = 0.4718 '- edgeBuffer
			highboundWv = 0.5032 '+ edgeBuffer
		ElseIf spectData.wvbandStr = "M4" Then
			'edgeBuffer = (0.571 - 0.5384)/(numWvPerBand*2)
			lowboundWv = 0.5384 '- edgeBuffer
			highboundWv = 0.571 '+ edgeBuffer
		ElseIf spectData.wvbandStr = "M5" Then
			'edgeBuffer = (0.687 - 0.654)/(numWvPerBand*2)
			lowboundWv = 0.654 '- edgeBuffer
			highboundWv = 0.687 '+ edgeBuffer
		ElseIf spectData.wvbandStr = "M6" Then
			'edgeBuffer = (0.7574 - 0.7322)/(numWvPerBand*2)
			lowboundWv = 0.7322 '- edgeBuffer
			highBoundWv = 0.7574 '+ edgeBuffer
		ElseIf spectData.wvbandStr = "M7" Then
			'edgeBuffer = (0.8892 - 0.8376)/(numWvPerBand*2)
			lowBoundWv = 0.8376 '- edgeBuffer
			highBoundWv = 0.8892 '+ edgeBuffer
		ElseIf spectData.wvbandStr = "I1" Then
			lowBoundWv =  0.5914 '- edgeBuffer
			highboundWv = 0.6874 '+ edgeBuffer
		ElseIf spectData.wvbandStr = "I2" Then
			lowboundWv = 0.8376 '- edgeBuffer
			highBoundWv = 0.8892 '+ edgeBuffer
		Else
			Print "Waveband not currently supported for definition with 'Spectra'."
		End If

		' Find the sources
		source0ID   = FindFullName("Optical Sources.Collimated Test Beam (0 deg pol)")
		source15ID  = FindFullName("Optical Sources.Collimated Test Beam (15 deg pol)")
		source30ID  = FindFullName("Optical Sources.Collimated Test Beam (30 deg pol)")
		source45ID  = FindFullName("Optical Sources.Collimated Test Beam (45 deg pol)")
		source60ID  = FindFullName("Optical Sources.Collimated Test Beam (60 deg pol)")
		source75ID  = FindFullName("Optical Sources.Collimated Test Beam (75 deg pol)")
		source90ID  = FindFullName("Optical Sources.Collimated Test Beam (90 deg pol)")
		source105ID = FindFullName("Optical Sources.Collimated Test Beam (105 deg pol)")
		source120ID = FindFullName("Optical Sources.Collimated Test Beam (120 deg pol)")
		source135ID = FindFullName("Optical Sources.Collimated Test Beam (135 deg pol)")
		source150ID = FindFullName("Optical Sources.Collimated Test Beam (150 deg pol)")
		source165ID = FindFullName("Optical Sources.Collimated Test Beam (165 deg pol)")
		sourceRID = FindFullName("Optical Sources.Collimated Test Beam (RCP)")
		sourceLID = FindFullName("Optical Sources.Collimated Test Beam (LCP)")

		' Make source wavelength definition with a spectrum
		SetSourceSpectrum source0ID, spectData.spectID
		SetSourceSpectrum source45ID, spectData.spectID
		SetSourceSpectrum source90ID, spectData.spectID
		SetSourceSpectrum source135ID, spectData.spectID
		SetSourceSpectrum sourceRID, spectData.spectID
		SetSourceSpectrum sourceLID, spectData.spectID

		' Options are that wavelengths 0: "As specified by list", 1: "Randomly according to spectrum", 2: "Evenly spaced, weighted according to spectrum"
		SetSourceRayWavelengthMethod source0ID, 2
		SetSourceRayWavelengthMethod source45ID, 2
		SetSourceRayWavelengthMethod source90ID, 2
		SetSourceRayWavelengthMethod source135ID, 2
		SetSourceRayWavelengthMethod sourceRID, 2
		SetSourceRayWavelengthMethod sourceLID, 2

		' Set wavelength range with upper and lower bounds
		SetSourceRayWavelengthRange source0ID, lowboundWv, highboundWv
		SetSourceRayWavelengthRange source45ID, lowboundWv, highboundWv
		SetSourceRayWavelengthRange source90ID, lowboundWv, highboundWv
		SetSourceRayWavelengthRange source135ID, lowboundWv, highboundWv
		SetSourceRayWavelengthRange sourceRID, lowboundWv, highboundWv
		SetSourceRayWavelengthRange sourceLID, lowboundWv, highboundWv

		' Set number of wavelengths (SetSourceWavelengthCount function doesn't work yet!)
		'SetSourceWavelengthCount source0ID, UBound(spectData.wavelengths)

	ElseIf spectData.spectTyp = List Then

		' More initialization
		' Get source(s)
		'sourceID = FindFullName("Optical Sources.Collimated Test Beam")
		source0ID   = FindFullName("Optical Sources.Collimated Test Beam (0 deg pol)")
		'source15ID  = FindFullName("Optical Sources.Collimated Test Beam (15 deg pol)")
		'source30ID  = FindFullName("Optical Sources.Collimated Test Beam (30 deg pol)")
		source45ID  = FindFullName("Optical Sources.Collimated Test Beam (45 deg pol)")
		'source60ID  = FindFullName("Optical Sources.Collimated Test Beam (60 deg pol)")
		'source75ID  = FindFullName("Optical Sources.Collimated Test Beam (75 deg pol)")
		source90ID  = FindFullName("Optical Sources.Collimated Test Beam (90 deg pol)")
		'source105ID = FindFullName("Optical Sources.Collimated Test Beam (105 deg pol)")
		'source120ID = FindFullName("Optical Sources.Collimated Test Beam (120 deg pol)")
		source135ID = FindFullName("Optical Sources.Collimated Test Beam (135 deg pol)")
		'source150ID = FindFullName("Optical Sources.Collimated Test Beam (150 deg pol)")
		'source165ID = FindFullName("Optical Sources.Collimated Test Beam (165 deg pol)")
		sourceRID = FindFullName("Optical Sources.Collimated Test Beam (RCP)")
		sourceLID = FindFullName("Optical Sources.Collimated Test Beam (LCP)")
		'sourceMMID = FindFullName("Optical Sources.Source for Mueller Matrix Sampling")

		' Options are that wavelengths 0: "As specified by list", 1: "Randomly according to spectrum", 2: "Evenly spaced, weighted according to spectrum"
		SetSourceRayWavelengthMethod source0ID, 0
		SetSourceRayWavelengthMethod source45ID, 0
		SetSourceRayWavelengthMethod source90ID, 0
		SetSourceRayWavelengthMethod source135ID, 0
		SetSourceRayWavelengthMethod sourceRID, 0
		SetSourceRayWavelengthMethod sourceLID, 0

		' Set their wavelengths
		' First clear out the existing wavelengths
		numWvForSource = GetSourceWavelengthCount(source0ID)
		For i = numWvForSource -1 To 1 Step -1
			DeleteSourceIthWavelengthSpec source0ID, i
			'DeleteSourceIthWavelengthSpec source15ID, i
			'DeleteSourceIthWavelengthSpec source30ID, i
			DeleteSourceIthWavelengthSpec source45ID, i
			'DeleteSourceIthWavelengthSpec source60ID, i
			'DeleteSourceIthWavelengthSpec source75ID, i
			DeleteSourceIthWavelengthSpec source90ID, i
			'DeleteSourceIthWavelengthSpec source105ID, i
			'DeleteSourceIthWavelengthSpec source120ID, i
			DeleteSourceIthWavelengthSpec source135ID, i
			'DeleteSourceIthWavelengthSpec source150ID, i
			'DeleteSourceIthWavelengthSpec source165ID, i
			DeleteSourceIthWavelengthSpec sourceRID, i
			DeleteSourceIthWavelengthSpec sourceLID, i
			'DeleteSourceIthWavelengthSpec sourceMMID, i
		Next i

		' Now set the new wavelengths
		wvSpec.weight = 1/UBound(spectData.wavelengths)
		For i = 0 To UBound(spectData.wavelengths) - 1
			If i = 0 Then
				wvSpec.wavelen = spectData.wavelengths(0)
				SetSourceIthWavelengthSpec source0ID, 0, wvSpec
				SetSourceIthWavelengthSpec source45ID, 0, wvSpec
				SetSourceIthWavelengthSpec source90ID, 0, wvSpec
				SetSourceIthWavelengthSpec source135ID, 0, wvSpec
				SetSourceIthWavelengthSpec sourceRID, 0, wvSpec
				SetSourceIthWavelengthSpec sourceLID, 0, wvSpec
				'SetSourceIthWavelengthSpec sourceMMID, 0, wvSpec
			Else
				wvSpec.wavelen = spectData.wavelengths(i)
				AddSourceWavelengthSpec source0ID, wvSpec
				AddSourceWavelengthSpec source45ID, wvSpec
				AddSourceWavelengthSpec source0ID, wvSpec
				AddSourceWavelengthSpec source0ID, wvSpec
				AddSourceWavelengthSpec source0ID, wvSpec
				AddSourceWavelengthSpec source0ID, wvSpec
				AddSourceWavelengthSpec source0ID, wvSpec
				'SetSourceIthWavelengthSpec sourceMMID, 0, wvSpec

			End If
		Next i
	Else
		Print "Invalid Wavelength definition!  Code currently supports 'list' and 'spectra'."

	End If

	' Set their directions
	SetSourceDirection source0ID, Sin(DegToRad(alpha)), Sin(DegToRad(beta))*Cos(DegToRad(alpha)), Cos(DegToRad(beta))*Cos(DegToRad(alpha))
	'SetSourceDirection source15ID, Sin(DegToRad(alpha)), Sin(DegToRad(beta))*Cos(DegToRad(alpha)), Cos(DegToRad(beta))*Cos(DegToRad(alpha))
	'SetSourceDirection source30ID, Sin(DegToRad(alpha)), Sin(DegToRad(beta))*Cos(DegToRad(alpha)), Cos(DegToRad(beta))*Cos(DegToRad(alpha))
	SetSourceDirection source45ID, Sin(DegToRad(alpha)), Sin(DegToRad(beta))*Cos(DegToRad(alpha)), Cos(DegToRad(beta))*Cos(DegToRad(alpha))
	'SetSourceDirection source60ID, Sin(DegToRad(alpha)), Sin(DegToRad(beta))*Cos(DegToRad(alpha)), Cos(DegToRad(beta))*Cos(DegToRad(alpha))
	'SetSourceDirection source75ID, Sin(DegToRad(alpha)), Sin(DegToRad(beta))*Cos(DegToRad(alpha)), Cos(DegToRad(beta))*Cos(DegToRad(alpha))
	SetSourceDirection source90ID, Sin(DegToRad(alpha)), Sin(DegToRad(beta))*Cos(DegToRad(alpha)), Cos(DegToRad(beta))*Cos(DegToRad(alpha))
	'SetSourceDirection source105ID, Sin(DegToRad(alpha)), Sin(DegToRad(beta))*Cos(DegToRad(alpha)), Cos(DegToRad(beta))*Cos(DegToRad(alpha))
	'SetSourceDirection source120ID, Sin(DegToRad(alpha)), Sin(DegToRad(beta))*Cos(DegToRad(alpha)), Cos(DegToRad(beta))*Cos(DegToRad(alpha))
	SetSourceDirection source135ID, Sin(DegToRad(alpha)), Sin(DegToRad(beta))*Cos(DegToRad(alpha)), Cos(DegToRad(beta))*Cos(DegToRad(alpha))
	'SetSourceDirection source150ID, Sin(DegToRad(alpha)), Sin(DegToRad(beta))*Cos(DegToRad(alpha)), Cos(DegToRad(beta))*Cos(DegToRad(alpha))
	'SetSourceDirection source165ID, Sin(DegToRad(alpha)), Sin(DegToRad(beta))*Cos(DegToRad(alpha)), Cos(DegToRad(beta))*Cos(DegToRad(alpha))
	'SetSourceDirection sourceMMID, Sin(DegToRad(alpha)), Sin(DegToRad(beta))*Cos(DegToRad(alpha)), Cos(DegToRad(beta))*Cos(DegToRad(alpha))
	If Not skipDOCP Then
		SetSourceDirection sourceRID, Sin(DegToRad(alpha)), Sin(DegToRad(beta))*Cos(DegToRad(alpha)), Cos(DegToRad(beta))*Cos(DegToRad(alpha))
		SetSourceDirection sourceLID, Sin(DegToRad(alpha)), Sin(DegToRad(beta))*Cos(DegToRad(alpha)), Cos(DegToRad(beta))*Cos(DegToRad(alpha))
	End If

	' Set their numbers of rays
	sourceDiameter = 7	' used on all analyses prior to 20140212
	'sourceDiameter = 8.4	' used for high-fidelity RSR sims, 20140212
	SetSourcePosGrid source0ID, sourceDiameter/2.0, sourceDiameter/2.0, numRays, numRays, False
	SetSourcePosGrid source45ID, sourceDiameter/2.0, sourceDiameter/2.0, numRays, numRays, False
	SetSourcePosGrid source90ID, sourceDiameter/2.0, sourceDiameter/2.0, numRays, numRays, False
	SetSourcePosGrid source135ID, sourceDiameter/2.0, sourceDiameter/2.0, numRays, numRays, False
	If Not skipDOCP Then
		SetSourcePosGrid sourceRID, sourceDiameter/2.0, sourceDiameter/2.0, numRays, numRays, False
		SetSourcePosGrid sourceLID, sourceDiameter/2.0, sourceDiameter/2.0, numRays, numRays, False
	End If
End Sub

Sub computePolarizationStatisticsTwelveSources(detID As Long, pathID As Long, DOLP As Double, modulation() As Double)

	' This script will compute the analyzer vector (i.e. the top row of the Mueller matrix) and the 
	' Degree of Linear Polarization (DOLP) at a detector by setting the polarization state of the input 
	' source, tracing rays, and determining the power on the detector. It assumes that twelve sources are 
	' traced, one for each polarization state. This is faster than changing the polarization state 
	' of a single source since it reduces the number of updates that must be made in FRED.

	Dim S(4) As Double, P(12) As Double, sourceAngles(4) As Double
	Dim diattenAngle As Double, alamda As Double
	Dim maxP As Double, minP As Double, chisq As Double, a(3) As Double, sig(4) As Double
	Dim source0ID As Long, source45ID As Long, source90ID As Long, source135ID As Long, i As Long
	Dim source15ID As Long, source30ID As Long, source60ID As Long, source75ID As Long, source105ID As Long
	Dim source120ID As Long, source150ID As Long, source165ID As Long
	Dim sourceRID As Long, sourceLID As Long, rayID As Long
	Dim ndata As Integer, ia(3) As Integer, ma As Integer
	Dim getRaySuccess As Boolean
	Dim rayData As T_RAY
	Dim eigenPolarization(4) As Double, DOP As Double, ellipticity As Double, ellipseAngle As Double
	Dim ellipseIsRightHanded As Boolean
	Dim adv As T_ADVANCEDRAYTRACE

	' Initialization	
	sourceAngles(0) = 0
	sourceAngles(1) = 45
	sourceAngles(2) = 90
	sourceAngles(3) = 135
	ndata = 4
	ma = 3
	alamda = -1
	ia(0) = 1	' Fit for all coefs
	ia(1) = 1
	ia(2) = 1
	eigenPolarization(0) = 1.
	InitAdvancedRaytrace adv
	adv.method = "User Path"	' use sequential raytracing to speed things up
	adv.pathID = pathID
	adv.draw = False
		
	' Get the source IDs
	source0ID   = FindFullName("Optical Sources.Collimated Test Beam (0 deg pol)")
	source15ID  = FindFullName("Optical Sources.Collimated Test Beam (15 deg pol)")
	source30ID  = FindFullName("Optical Sources.Collimated Test Beam (30 deg pol)")
	source45ID  = FindFullName("Optical Sources.Collimated Test Beam (45 deg pol)")
	source60ID  = FindFullName("Optical Sources.Collimated Test Beam (60 deg pol)")
	source75ID  = FindFullName("Optical Sources.Collimated Test Beam (75 deg pol)")
	source90ID  = FindFullName("Optical Sources.Collimated Test Beam (90 deg pol)")
	source105ID = FindFullName("Optical Sources.Collimated Test Beam (105 deg pol)")
	source120ID = FindFullName("Optical Sources.Collimated Test Beam (120 deg pol)")
	source135ID = FindFullName("Optical Sources.Collimated Test Beam (135 deg pol)")
	source150ID = FindFullName("Optical Sources.Collimated Test Beam (150 deg pol)")
	source165ID = FindFullName("Optical Sources.Collimated Test Beam (165 deg pol)")
	
	' Do the raytrace
	EnableTextPrinting False
	AdvancedRaytrace adv	
	EnableTextPrinting True

	' Loop over the rays on the detector and accumulate the power from each source
	P(0) = 0	' 0 degrees
	P(1) = 0	' 15 degrees
	P(2) = 0	' 30 degrees
	P(3) = 0	' 45 degrees
	P(4) = 0	' 60 degrees
	P(5) = 0	' 75 degrees
	P(6) = 0	' 90 degrees
	P(7) = 0	' 105 degrees
	P(8) = 0	' 120 degrees
	P(9) = 0	' 135 degrees
	P(10) = 0	' 150 degrees
	P(11) = 0	' 165 degrees
	getRaySuccess = GetFirstRay(rayID, rayData)
	While getRaySuccess
		If rayData.entity = detID Then
			' Which source did this ray come from?
			If rayData.Source = source0ID Then
				P(0) = P(0) + rayData.power
			ElseIf rayData.Source = source15ID Then
				P(1) = P(1) + rayData.power
			ElseIf rayData.Source = source30ID Then
				P(2) = P(2) + rayData.power
			ElseIf rayData.Source = source45ID Then
				P(3) = P(3) + rayData.power
			ElseIf rayData.Source = source60ID Then
				P(4) = P(4) + rayData.power
			ElseIf rayData.Source = source75ID Then
				P(5) = P(5) + rayData.power
			ElseIf rayData.Source = source90ID Then
				P(6) = P(6) + rayData.power
			ElseIf rayData.Source = source105ID Then
				P(7) = P(7) + rayData.power
			ElseIf rayData.Source = source120ID Then
				P(8) = P(8) + rayData.power
			ElseIf rayData.Source = source135ID Then
				P(9) = P(9) + rayData.power
			ElseIf rayData.Source = source150ID Then
				P(10) = P(10) + rayData.power
			ElseIf rayData.Source = source165ID Then
				P(11) = P(11) + rayData.power
			End If
		End If			
		getRaySuccess = GetNextRay(rayID, rayData)
	Wend
	
	Print GetName(detID); "	"; P(0); "	";  P(1); "	";  P(2); "	";  P(3); "	";  P(4); "	";  P(5); "	";  P(6); "	";  P(7); "	";  P(8); "	";  P(9); "	";  P(10); "	";  P(11)
	
	'Print "P0 = ", P0, " Watts"
	'Print "P45 = ", P45, " Watts"
	'Print "P90 = ", P90, " Watts"
	'Print "P135 = ", P135, " Watts"
	'Print "PR = ", PR, " Watts"
	'Print "PL = ", PL, " Watts"
End Sub