!***************************************************************************
!
!								Index vs Wavelength
!
!***************************************************************************

SUBROUTINE INFO_Index_Wavelength
	USE PM_Data

	PlotName = 'Index of Refraction (Wavelength)'
	XName = 'Wavelength (µm)'
	YName = 'Nx;Ny;Nz'
	PlotSize(1) = 50
	PlotSize(2) = 3
	PlotSize(3) = 1
	Plot_Type = PLOT_2D
	XAxisIsAngle = .FALSE.
	YAxisIsAngle = .FALSE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Temperature_Pos) = 1
END SUBROUTINE INFO_Index_Wavelength

SUBROUTINE Plot_Index_Wavelength(XAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)
	
!	Local variables
	INTEGER XIndex
	REAL*8 Step, Lambda, Nx, Ny, Nz

	IF (AutoRange) THEN
		XMin = Abs_Min
		XMax = Abs_Max
	END IF

!	Set the step size
    Step = (XMax - XMin) / (XRes - 1)

!	Perform Calculations
	DO XIndex = 1, XRes, 1
		Lambda = XMin + (XIndex - 1)*Step 
		CALL INDEX(Crystal, Lambda, Temperature,Nx, Ny, Nz)
		XAxis(XIndex) = Lambda
		PlotData(XIndex,1,1) = Nx
		PlotData(XIndex,2,1) = Ny
		PlotData(XIndex,3,1) = Nz
	ENDDO

END SUBROUTINE Plot_Index_Wavelength

!***************************************************************************
!
!								Index vs Temperature
!
!***************************************************************************

SUBROUTINE INFO_INDEX_TEMPERATURE
	USE PM_Data

	PlotName = 'Index of Refraction (Temperature)'
	XName = 'Temperature (C)'
	YName = 'Nx;Ny;Nz'
	PlotSize(1) = 50
	PlotSize(2) = 3
	PlotSize(3) = 1
	Plot_Type = PLOT_2D
	XAxisIsAngle = .FALSE.
	YAxisIsAngle = .FALSE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
END SUBROUTINE INFO_INDEX_TEMPERATURE

SUBROUTINE Plot_INDEX_TEMPERATURE(XAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)
	
!	Local variables
	INTEGER XIndex
	REAL*8 Step, Lambda, Nx, Ny, Nz

	IF (AutoRange) THEN
		XMin = 20
		XMax = 180
	END IF

!	Set the step size
    Step = (XMax - XMin) / (XRes - 1)

!	Perform Calculations
	DO XIndex = 1, XRes, 1
		Temperature = XMin + (XIndex - 1)*Step 
		CALL INDEX(Crystal, Lambda_Pump, Temperature,Nx, Ny, Nz)
		XAxis(XIndex) = Temperature
		PlotData(XIndex,1,1) = Nx
		PlotData(XIndex,2,1) = Ny
		PlotData(XIndex,3,1) = Nz
	ENDDO

END SUBROUTINE Plot_INDEX_TEMPERATURE

!***************************************************************************
!
!              Slow Index - Fast Index (Phi, Theta)
!
!***************************************************************************

SUBROUTINE INFO_DIFFN_PHI_THETA
	USE PM_Data

	PlotName = 'Slow Index - Fast Index (Phi, Theta)'
	XName = 'Phi Pump (deg)'
	YName = 'Theta Pump (deg)'
	PlotSize(1) = 71
	PlotSize(2) = 51
	PlotSize(3) = 1
	Plot_Type = PLOT_3D
	XAxisIsAngle = .TRUE.
	YAxisIsAngle = .TRUE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
END SUBROUTINE INFO_DIFFN_PHI_THETA

SUBROUTINE PLOT_DIFFN_PHI_THETA(XAxis, YAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)
	REAL*8, INTENT(INOUT) :: YAxis(YRes)
	
!	Local variables
	INTEGER  XIndex, YIndex	
	REAL*8 NF, NS, Min_Lim, Max_Lim
	REAL*8 Theta_Step, Phi_Step

!	Set range
	IF (AutoRange) THEN
		XMin = 0d0
		XMax = 2d0 * Pi
		YMin = 0d0
		YMax = Pi
	END IF

	Phi_Step = (XMax - XMin) / (XRes - 1)
	Theta_Step = (YMax - YMin) / (YRes - 1)

!	Perform Calculations
	DO XIndex = 1, XRes, 1
		Phi_Pump = XMin + (XIndex - 1)*Phi_Step 
			
		DO YIndex = 1, YRes, 1
			Theta_Pump = YMin + (YIndex - 1)*Theta_Step
	
			CALL GET_INDICES(NF,NS,Crystal,Lambda_Pump,Temperature,Theta_Pump,Phi_Pump)

			XAxis(XIndex) = Phi_Pump * 180 / Pi
			YAxis(YIndex) = Theta_Pump * 180 / Pi
			PlotData(XIndex,YIndex,1) = NS - NF
		ENDDO
	ENDDO	
	IF ((Lambda_Pump .GT. Abs_Max) .OR. (Lambda_Pump .LT. Abs_Min)) THEN
		WRITE(*,*) 'WARNING: The pump wavelength is not valid for this crystal.  Results may be erroneous'
		StatusOK = .FALSE.
	ENDIF
END SUBROUTINE PLOT_DIFFN_PHI_THETA

!***************************************************************************
!
!              Maximum Phase Match Function (Phi and Theta)
!
!***************************************************************************

SUBROUTINE INFO_MAXPMF_PHI_THETA
	USE PM_Data

	PlotName = 'Maximum Phase Match Function (Phi and Theta)'
	XName = 'Phi Pump (deg)'
	YName = 'Theta Pump (deg)'
	PlotSize(1) = 71
	PlotSize(2) = 51
	PlotSize(3) = 1
	Plot_Type = PLOT_3D
	XAxisIsAngle = .TRUE.
	YAxisIsAngle = .TRUE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Pump_Width_Pos) = 1
	Params_Used(Crystal_Length_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Lambda_Sgnl_Pos) = 1
	Params_Used(TypePM_Pos)      = 1
	Params_Used(Use_Previous_Pos)= 1
END SUBROUTINE INFO_MAXPMF_PHI_THETA

SUBROUTINE PLOT_MAXPMF_PHI_THETA(XAxis, YAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)
	REAL*8, INTENT(INOUT) :: YAxis(YRes)

!	Local variables
	INTEGER Info, XIndex, YIndex
	REAL*8 Theta_Step, Phi_Step, PMF_Value, Guess(2)
	LOGICAL WAVELENGTHS_VALID, Success
		
!	Set range
	IF (AutoRange) THEN
		XMin = 0d0
		XMax = 2d0 * Pi
		YMin = 0d0
		YMax = Pi
	END IF
	Phi_Step = (XMax - XMin) / (XRes - 1)
	Theta_Step = (YMax - YMin) / (YRes - 1)

!	Perform Calculations
	Guess(1)=Theta_Sgnl_Guess
	Guess(2)=Theta_Idlr_Guess
	DO XIndex = 1, XRes, 1
		Phi_Pump = XMin + (XIndex - 1)*Phi_Step 
		WRITE(*,"(I3,'% Done')") FLOOR(100.0d0*XIndex/XRes)
	
		DO YIndex = 1, YRes, 1
			Theta_Pump = YMin + (YIndex - 1)*Theta_Step

			CALL INIT_NPUMP 

			CALL FIND_OPT_PMF_2D(Guess, PMF_Value, Info)
						
			IF ( Use_Previous ) THEN
				Guess(1) = Theta_Sgnl
				Guess(2) = Theta_Idlr
			ELSE
				Guess(1) = Theta_Sgnl_Guess
				Guess(2) = Theta_Idlr_Guess
			ENDIF
			PlotData(XIndex,YIndex,1) = PMF_Value
			XAxis(XIndex) = Phi_Pump * 180d0 / Pi
			YAxis(YIndex) = Theta_Pump * 180d0 / Pi
		ENDDO
	ENDDO
	IF (.NOT. WAVELENGTHS_VALID() ) THEN
		WRITE(*,*) 'WARNING: Wavelength out of range.  Results may not be valid' 
		WRITE(*,"('Pump:',F8.4,'; Signal:',F8.4,'; Idler:',F8.4)") Lambda_Pump, Lambda_Sgnl, Lambda_Idlr
		StatusOK = .FALSE.
	ENDIF
END SUBROUTINE PLOT_MAXPMF_PHI_THETA

!***************************************************************************
!
!     Phase Match Function For Colinear Downconversion (Theta_Pump,Phi_Pump)
!
!***************************************************************************

SUBROUTINE INFO_PMFCOL_PHI_THETA
	USE PM_Data

	PlotName = 'Phase Match Function For Colinear Downconversion (Phi_Pump,Theta_Pump)'
	XName = 'Phi Pump (deg)'
	YName = 'Theta Pump (deg)'
	PlotSize(1) = 500
	PlotSize(2) = 500
	PlotSize(3) = 1
	Plot_Type = PLOT_3D
	XAxisIsAngle = .TRUE.
	YAxisIsAngle = .TRUE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Pump_Width_Pos) = 1
	Params_Used(Crystal_Length_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Lambda_Sgnl_Pos) = 1
	Params_Used(Phi_Sgnl_Pos) = 1
	Params_Used(TypePM_Pos) = 1
	Params_Used(Use_Previous_Pos) = 1
END SUBROUTINE INFO_PMFCOL_PHI_THETA

SUBROUTINE PLOT_PMFCOL_PHI_THETA(XAxis, YAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)
	REAL*8, INTENT(INOUT) :: YAxis(YRes)

!	Local variables
	INTEGER Info, XIndex, YIndex
	REAL*8 PMF_Value, Phi_Step, Theta_Step, PHASE_MATCH_FUNC
	LOGICAL WAVELENGTHS_VALID, Success
	        
!	Set range
	IF (AutoRange) THEN
		XMin = 0d0
		XMax = 2d0 * Pi
		YMin = 0d0
		YMax = Pi
	END IF

	Phi_Step   = (XMax - XMin) / (XRes - 1)
	Theta_Step = (YMax - YMin) / (YRes - 1)

	Theta_Sgnl = 0d0
	Theta_Idlr = 0d0
	
	Wavelength_Violation = .FALSE.
!	Perform calculations
	DO XIndex = 1, XRes, 1
		WRITE(*,"(I3,'% Done')") FLOOR(100.0d0*XIndex/XRes)
		Phi_Pump = XMin + (XIndex - 1)*Phi_Step 

		DO YIndex = 1, YRes, 1
			Theta_Pump = YMin + (YIndex - 1) * Theta_Step
			CALL INIT_NPUMP

			XAxis(XIndex) = Phi_Pump * 180d0 / Pi
			YAxis(YIndex) = Theta_Pump * 180d0 / Pi
			PlotData(XIndex,YIndex,1) = PHASE_MATCH_FUNC(Success)
		ENDDO

		IF ( .NOT. WAVELENGTHS_VALID() ) THEN
			WRITE(*,"('WAVELENGTH VIOLATION - Pump:',F8.4,'; Signal:',F8.4,'; Idler:',F8.4)") Lambda_Pump, Lambda_Sgnl, Lambda_Idlr
			Wavelength_Violation = .TRUE.
		ENDIF
	ENDDO

END SUBROUTINE PLOT_PMFCOL_PHI_THETA

!***************************************************************************
!
!     Phase Match Function For Colinear Downconversion (Lambda,Theta_Pump)
!
!***************************************************************************

SUBROUTINE INFO_PMFCOL_LAMBDA_THETA
	USE PM_Data

	PlotName = 'Phase Match Function For Colinear Downconversion (Lambda,Theta_Pump)'
	XName = 'Lambda Signal (µm)'
	YName = 'Theta Pump (deg)'
	PlotSize(1) = 500
	PlotSize(2) = 500
	PlotSize(3) = 1
	Plot_Type = PLOT_3D
	XAxisIsAngle = .FALSE.
	YAxisIsAngle = .TRUE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Pump_Width_Pos) = 1
	Params_Used(Crystal_Length_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Phi_Pump_Pos) = 1
	Params_Used(Phi_Sgnl_Pos) = 1
	Params_Used(TypePM_Pos) = 1
	Params_Used(Use_Previous_Pos) = 1
END SUBROUTINE INFO_PMFCOL_LAMBDA_THETA

SUBROUTINE PLOT_PMFCOL_LAMBDA_THETA(XAxis, YAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)
	REAL*8, INTENT(INOUT) :: YAxis(YRes)

!	Local variables
	INTEGER Info, XIndex, YIndex
	REAL*8 PMF_Value, Lambda_Step, Theta_Step, PHASE_MATCH_FUNC
	LOGICAL WAVELENGTHS_VALID, Success
	        
!	Set range
	IF (AutoRange) THEN
		XMin = 1/(1/Lambda_Pump - 1/Abs_Max) 	
		IF (XMin .LT. Abs_Min) THEN 
			XMin = Abs_Min
		END IF
		XMax = Abs_Max
!		It doesn't like the wavelength right on the edge, so move in a bit
		XMin = XMin+(XMax-XMin)/XRes
		YMin = 0d0
		YMax = Pi
	END IF

	Lambda_Step = (XMax - XMin) / (XRes - 1)
	Theta_Step  = (YMax - YMin) / (YRes - 1)

	Theta_Sgnl = 0d0
	Theta_Idlr = 0d0
	Wavelength_Violation = .FALSE.
!	Perform calculations
	DO XIndex = 1, XRes, 1
		WRITE(*,"(I3,'% Done')") FLOOR(100.0d0*XIndex/XRes)
		Lambda_Sgnl = XMin + (XIndex - 1)*Lambda_Step 
		Lambda_Idlr = 1/(1/Lambda_Pump - 1/Lambda_Sgnl)

		DO YIndex = 1, YRes, 1
			Theta_Pump = YMin + (YIndex - 1) * Theta_Step
			CALL INIT_NPUMP

			XAxis(XIndex) = Lambda_Sgnl
			YAxis(YIndex) = Theta_Pump * 180d0 / Pi
			PlotData(XIndex,YIndex,1) = PHASE_MATCH_FUNC(Success)
		ENDDO

		IF ( .NOT. WAVELENGTHS_VALID() ) THEN
			WRITE(*,"('WAVELENGTH VIOLATION - Pump:',F8.4,'; Signal:',F8.4,'; Idler:',F8.4)") Lambda_Pump, Lambda_Sgnl, Lambda_Idlr
			Wavelength_Violation = .TRUE.
		ENDIF
	ENDDO

END SUBROUTINE PLOT_PMFCOL_LAMBDA_THETA


!***************************************************************************
!
!              Optimum Theta Signal at Chosen Phi (Lambda)
!
!***************************************************************************

SUBROUTINE INFO_THETAOPT_LAMBDA
	USE PM_Data

	PlotName = 'Optimum Theta Signal at Chosen Phi (Lambda)'
	XName = 'Lambda Signal (µm)'
	YName = 'Theta Signal Internal (deg);Theta Signal External (deg);Theta of Corellated Photon (deg);Theta External of Correlated Photon(deg);D Effective'
	PlotSize(1) = 50
	PlotSize(2) = 5
	PlotSize(3) = 1
	Plot_Type = PLOT_2D
	XAxisIsAngle = .FALSE.
	YAxisIsAngle = .FALSE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Pump_Width_Pos) = 1
	Params_Used(Crystal_Length_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Theta_Pump_Pos) = 1
	Params_Used(Phi_Pump_Pos) = 1
	Params_Used(Phi_Sgnl_Pos) = 1
	Params_Used(TypePM_Pos) = 1
	Params_Used(Use_Previous_Pos) = 1
END SUBROUTINE INFO_THETAOPT_LAMBDA

SUBROUTINE PLOT_THETAOPT_LAMBDA(XAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)

!	Local variables
	INTEGER Info, XIndex
	REAL*8 PMF_Value, Guess(2), NF,NS, Lambda_Step, GET_D_EFF
	REAL*8 Nx, Ny, Nz, Theta_Sgnl_Ext, Theta_Idlr_Ext, Argument
	LOGICAL WAVELENGTHS_VALID
	        
!	Set range
	IF (AutoRange) THEN
		XMin = 1/(1/Lambda_Pump - 1/Abs_Max) 	
		IF (XMin .LT. Abs_Min) THEN 
			XMin = Abs_Min
		END IF
		XMax = Abs_Max
!		It doesn't like the wavelength right on the edge, so move in a bit
		XMin = XMin+(XMax-XMin)/XRes
	END IF

	Lambda_Step = (XMax - XMin) / (XRes - 1)

	Theta_Sgnl_Guess = 0
	Theta_Idlr_Guess = 0

!	Perform calculations
	Guess(1) = Theta_Sgnl_Guess
	Guess(2) = Theta_Idlr_Guess 
	Wavelength_Violation = .FALSE.
	DO XIndex = 1, XRes, 1
		WRITE(*,"(I3,'% Done')") FLOOR(100.0d0*XIndex/XRes)
		Lambda_Sgnl = XMin + (XIndex - 1)*Lambda_Step 
		Lambda_Idlr = 1/(1/Lambda_Pump - 1/Lambda_Sgnl)

!		Calculate theta_signal optimum and theta_idler optimum
!		IF ( TypePM .EQ. 1) THEN
!			CALL INDEX(Crystal, Lambda_Sgnl, Temperature, N_Sgnl, Ny, Nz)
!			CALL INDEX(Crystal, Lambda_Idlr, Temperature, N_Idlr, Ny, Nz)
!
!			Argument = 0.5d0*Lambda_Pump / (N_Pump*N_Sgnl*Lambda_Sgnl) *	&
!						( N_Sgnl**2 + (N_Pump*Lambda_Sgnl/Lambda_Pump)**2	&
!						 - (N_Idlr*Lambda_Sgnl/Lambda_Idlr)**2 )			
!			IF (DABS(Argument) .LE. 1) THEN
!				Theta_Sgnl = DACOS(Argument)
!			ELSE
!				Theta_Sgnl = 0
!			ENDIF
!							
!			Argument = 0.5d0*Lambda_Pump / (N_Pump*N_Idlr*Lambda_Idlr) *	&
!						( N_Idlr**2 + (N_Pump*Lambda_Idlr/Lambda_Pump)**2	&
!						 - (N_Sgnl*Lambda_Idlr/Lambda_Sgnl)**2 )			
!			IF (DABS(Argument) .LE. 1) THEN
!				Theta_Idlr = DACOS(Argument)
!			ELSE
!				Theta_Idlr = 0
!			ENDIF
!		ELSE
		CALL FIND_OPT_PMF_2D(Guess,PMF_Value,Info)

!		Give preference to keeping the answer the same sign as the guess
		IF ( DSIGN(1d0,Theta_Sgnl) .NE. DSIGN(1d0,Theta_Sgnl_Guess) ) THEN
			Guess(1) = DSIGN(Theta_Sgnl,Guess(1))
			Guess(2) = DSIGN(Theta_Idlr,Guess(2))
			CALL FIND_OPT_PMF_2D(Guess,PMF_Value,Info)
		ENDIF

!		Initial guesses for theta (signal,idler) optimum in next iteration
		IF (Use_Previous .AND. (PMF_Value .GT. 0.05d0)) THEN
			Guess(1) = Theta_Sgnl
			Guess(2) = Theta_Idlr
		ELSE
			Guess(1) = Theta_Sgnl_Guess
			Guess(2) = Theta_Idlr_Guess
		ENDIF

		Argument = DSIN( Theta_Sgnl )*N_Sgnl
		IF (DABS(Argument) .LE. 1.0d0) THEN
			Theta_Sgnl_Ext = DASIN( Argument ) 
		ELSE
			IF (Theta_Sgnl .GT. 0d0) THEN
				Theta_Sgnl_Ext = Pi-Theta_Sgnl
			ELSE
				Theta_Sgnl_Ext = -Pi+Theta_Sgnl
			ENDIF						
		ENDIF
		Argument = DSIN( Theta_Idlr )*N_Idlr
		IF (DABS(Argument) .LE. 1.0d0) THEN
			Theta_Idlr_Ext = DASIN( Argument ) 
		ELSE
			IF (Theta_Idlr .GT. 0d0) THEN
				Theta_Idlr_Ext = Pi-Theta_Idlr
			ELSE
				Theta_Idlr_Ext = -Pi+Theta_Idlr
			ENDIF						
		ENDIF

		XAxis(XIndex) = Lambda_Sgnl
		PlotData(XIndex,1,1) = Theta_Sgnl * 180/Pi
		PlotData(XIndex,2,1) = Theta_Sgnl_Ext * 180/Pi
		PlotData(XIndex,3,1) = Theta_Idlr *180/Pi
		PlotData(XIndex,4,1) = Theta_Idlr_Ext *180d0/Pi
		PlotData(XIndex,5,1) = GET_D_EFF()
		IF ( .NOT. WAVELENGTHS_VALID() ) THEN
			WRITE(*,"('WAVELENGTH VIOLATION - Pump:',F8.4,'; Signal:',F8.4,'; Idler:',F8.4)") Lambda_Pump, Lambda_Sgnl, Lambda_Idlr
			Wavelength_Violation = .TRUE.
		ENDIF
	ENDDO

END SUBROUTINE PLOT_THETAOPT_LAMBDA

!***************************************************************************
!
!              Phase Matching Function (Lambda, Theta)
!
!***************************************************************************

SUBROUTINE INFO_PMF_LAMBDA_THETA
	USE PM_Data

	PlotName = 'Phase Matching Function (Lambda, Theta)'
	XName = 'Lambda Signal (µm)'
	IF (Angle_External) THEN
		YName = 'Theta Signal External(deg)'
	ELSE
		YName = 'Theta Signal Internal(deg)'
	ENDIF
	PlotSize(1) = 50
	PlotSize(2) = 50
	PlotSize(3) = 2
	Plot_Type = PLOT_3D
	XAxisIsAngle = .FALSE.
	YAxisIsAngle = .TRUE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Pump_Width_Pos) = 1
	Params_Used(Crystal_Length_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Theta_Pump_Pos) = 1
	Params_Used(Phi_Pump_Pos) = 1
	Params_Used(Angle_External_Pos) = 1
	Params_Used(TypePM_Pos) = 1
END SUBROUTINE INFO_PMF_LAMBDA_THETA
 

SUBROUTINE PLOT_PMF_LAMBDA_THETA(XAxis, YAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)
	REAL*8, INTENT(INOUT) :: YAxis(YRes)

!	Local variables
	INTEGER Info, XIndex, YIndex
	LOGICAL WAVELENGTHS_VALID, Success
	REAL*8 Guess, Theta_Sgnl_Ext, PMF_Value, Argument
	REAL*8 Lambda_Step, Theta_Step, DEff

!	Set x-axis range
	IF (AutoRange) THEN
		XMin = 1/(1/Lambda_Pump - 1/Abs_Max) 
		XMax = Abs_Max
!		It doesn't like the wavelength right on the edge, so move in a bit
		XMin = XMin+(XMax-XMin)/XRes
		YMin = 0.0d0
		YMax = Pi/18d0
		YRes = 50
	END IF
	Lambda_Step = (XMax - XMin) / (XRes-1)
	Theta_Step = (YMax - YMin) / (YRes-1)
	
!	Perform calculations
	Guess = Theta_Idlr_Guess
	Wavelength_Violation = .FALSE.
	DO XIndex = 1, XRes, 1
		WRITE(*,"(I3,'% Done')") FLOOR(100.0d0*XIndex/XRes)

		Lambda_Sgnl = XMin + (XIndex - 1)*Lambda_Step 
		Lambda_Idlr=1/(1/Lambda_Pump-1/Lambda_Sgnl)

		DO YIndex = 1, YRes, 1
			IF (Angle_External) THEN
				Theta_Sgnl_Ext = YMin+(YIndex - 1)*Theta_Step
				CALL GetThetaInternal(Theta_Sgnl_Ext,Theta_Sgnl)
			ELSE
				Theta_Sgnl = YMin + (YIndex - 1)*Theta_Step
			ENDIF
				
			IF (YIndex .EQ. 1) THEN
				Guess = Theta_Sgnl
			ELSE
				Guess = Theta_Idlr
			ENDIF
					
			CALL FIND_OPT_PMF_1D(Guess, PMF_Value, Info)
					
			Argument = N_Sgnl * SIN(Theta_Sgnl)
			IF (DABS(Argument) .LE. 1.0d0) THEN
				Theta_Sgnl_Ext = DASIN( Argument ) 
			ELSE
				IF (Theta_Sgnl .GT. 0d0) THEN
					Theta_Sgnl_Ext = Pi-Theta_Sgnl
				ELSE
					Theta_Sgnl_Ext = -Pi+Theta_Sgnl
				ENDIF						
			ENDIF

			CALL Calculate_DEff(DEff)
			PlotData(XIndex,YIndex,1) = PMF_Value
			PlotData(XIndex,YIndex,2) = DEff
			XAxis(XIndex) = Lambda_Sgnl
			IF (Angle_External) THEN
				YAxis(YIndex) = Theta_Sgnl_Ext * 180.0d0/Pi
			ELSE
				YAxis(YIndex) = Theta_Sgnl * 180.0d0/Pi
			ENDIF
		ENDDO
		IF ( .NOT. WAVELENGTHS_VALID() ) THEN
			WRITE(*,"('WAVELENGTH VIOLATION - Pump:',F8.4,'; Signal:',F8.4,'; Idler:',F8.4)") Lambda_Pump, Lambda_Sgnl, Lambda_Idlr
			Wavelength_Violation = .TRUE.
		ENDIF
	ENDDO

END SUBROUTINE PLOT_PMF_LAMBDA_THETA

!***************************************************************************
!
!		SUBROUTINE PLOT_PMFINT_LAMBDA_THETA(PMVars, PlotInfo, PlotData, XAxis, YAxis)
!
!		Creates data for a 3D plot, Phase-matching function = f(Theta,Lambda)
!
!***************************************************************************

SUBROUTINE INFO_PMFINT_LAMBDA_THETA
	USE PM_DATA
	IMPLICIT NONE

	PlotName = 'Phase Matching Function Integral (Lambda, Theta)'
	XName = 'Lambda Signal (µm)'
	IF (Angle_External) THEN
		YName = 'Theta Signal External(deg)'
	ELSE
		YName = 'Theta Signal Internal(deg)'
	ENDIF
	PlotSize(1) = 50
	PlotSize(2) = 50
	PlotSize(3) = 1
	Plot_Type = PLOT_3D
	XAxisIsAngle = .FALSE.
	YAxisIsAngle = .TRUE.
	Params_Used(Pump_Width_Pos) = 1
	Params_Used(Crystal_Length_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Theta_Pump_Pos) = 1
	Params_Used(Phi_Pump_Pos) = 1
	Params_Used(Angle_External_Pos) = 1
	Params_Used(TypePM_Pos) = 1
END SUBROUTINE INFO_PMFINT_LAMBDA_THETA

SUBROUTINE PLOT_PMFINT_LAMBDA_THETA(XAxis, YAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)
	REAL*8, INTENT(INOUT) :: YAxis(YRes)

!	Local variables
	INTEGER Info, XIndex, YIndex
	REAL*8 Guess
	REAL*8 Theta_Sgnl_Ext, PMFInt_Value, Argument
	REAL*8 Nx, Ny, Nz, PMF_INTEGRAL
	LOGICAL WAVELENGTHS_VALID, Success
	REAL*8 Lambda_Step, Theta_Step, Max

!	Set x-axis range
	IF (AutoRange) THEN
		XMin = 1/(1/Lambda_Pump - 1/Abs_Max) 
		XMax = Abs_Max
!		It doesn't like the wavelength right on the edge, so move in a bit
		XMin = XMin+(XMax-XMin) / XRes
		YMin = 0.0d0
		YMax = Pi / 18d0
	ENDIF
	Lambda_Step = (XMax - XMin) / (XRes - 1)
	Theta_Step  = (YMax - YMin) / (YRes - 1)

!	Perform calculations
	Guess = Theta_Idlr_Guess
	DO XIndex = 1, XRes, 1
		WRITE(*,"(I3,'% Done')") FLOOR(100.0d0*XIndex/XRes)

		Lambda_Sgnl = XMin + (XIndex - 1)*Lambda_Step 
		Lambda_Idlr=1/(1/Lambda_Pump-1/Lambda_Sgnl)

		DO YIndex = 1, YRes, 1
			IF (Angle_External) THEN
				Theta_Sgnl_Ext = YMin+(YIndex - 1)*Theta_Step
				CALL GetThetaInternal(Theta_Sgnl_Ext,Theta_Sgnl)
			ELSE
				Theta_Sgnl = YMin + (YIndex - 1)*Theta_Step
			ENDIF
					
			IF (YIndex .EQ. 1) THEN
				Guess = Theta_Sgnl
			ELSE
				Guess = Theta_Idlr
			ENDIF
					
			PMFInt_Value = PMF_INTEGRAL(Success)
					
			Argument = N_Sgnl * SIN(Theta_Sgnl)
			IF (DABS(Argument) .LE. 1.0d0) THEN
				Theta_Sgnl_Ext = DASIN( Argument ) 
			ELSE
				IF (Theta_Sgnl .GT. 0d0) THEN
					Theta_Sgnl_Ext = Pi-Theta_Sgnl
				ELSE
					Theta_Sgnl_Ext = -Pi+Theta_Sgnl
				ENDIF						
			ENDIF

			PlotData(XIndex,YIndex,1) = PMFInt_Value
			XAxis(XIndex) = Lambda_Sgnl
			IF (Angle_External) THEN
				YAxis(YIndex) = Theta_Sgnl_Ext * 180.0d0/Pi
			ELSE
				YAxis(YIndex) = Theta_Sgnl * 180.0d0/Pi
			ENDIF
		ENDDO      
		IF ( .NOT. WAVELENGTHS_VALID() ) THEN
			WRITE(*,"('WAVELENGTH VIOLATION - Pump:',F8.4,'; Signal:',F8.4,'; Idler:',F8.4)") Lambda_Pump, Lambda_Sgnl, Lambda_Idlr
			Wavelength_Violation = .TRUE.
		ENDIF
	ENDDO

	DO XIndex = 1, XRes, 1
		Max = 0d0
		DO YIndex = 1, YRes, 1
			IF (PlotData(XIndex,YIndex,1) .GT. Max) THEN
				Max = PlotData(XIndex,YIndex,1)
			ENDIF
		ENDDO
		DO YIndex = 1, YRes, 1
			PlotData(XIndex,YIndex,1) = PlotData(XIndex,YIndex,1) / Max
		ENDDO
	ENDDO

END SUBROUTINE PLOT_PMFINT_LAMBDA_THETA

!***************************************************************************
!
!             Phase Matching Function (Lambda, Phi)
!
!***************************************************************************

SUBROUTINE INFO_PMF_LAMBDA_PHI
	USE PM_Data

	PlotName = 'Phase Matching Function (Lambda, Phi)'
	XName = 'Lambda Signal (µm)'
	YName = 'Phi Signal (deg)'
	PlotSize(1) = 50
	PlotSize(2) = 50
	PlotSize(3) = 1
	Plot_Type = PLOT_3D
	XAxisIsAngle = .FALSE.
	YAxisIsAngle = .TRUE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Pump_Width_Pos) = 1
	Params_Used(Crystal_Length_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Theta_Pump_Pos) = 1
	Params_Used(Phi_Pump_Pos) = 1
	Params_Used(TypePM_Pos) = 1
	Params_Used(Use_Previous_Pos) = 1
	Params_Used(Phi_Idlr_Pos) = 1
END SUBROUTINE INFO_PMF_LAMBDA_PHI

SUBROUTINE PLOT_PMF_LAMBDA_PHI(XAxis, YAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)
	REAL*8, INTENT(INOUT) :: YAxis(YRes)

!	Local variables
	INTEGER XIndex, YIndex, Info
	REAL*8 PMF_Value, PHASE_MATCH_FUNC, Lambda_Step, Phi_Step, Guess(2)
	LOGICAL WAVELENGTHS_VALID, Success

!	Set x-axis range
	IF (AutoRange) THEN
		XMin = 1/(1/Lambda_Pump - 1/Abs_Max) 
		XMax = Abs_Max
!		It doesn't like the wavelength right on the edge, so move in a bit
		XMin = XMin+(XMax-XMin)/XRes
		YMin = -0.5d0 * Pi/180d0
		YMax = 0.5d0 * Pi/180d0
		YRes = 50
	END IF
	Lambda_Step = (XMax - XMin) / (XRes - 1)
	Phi_Step = (YMax - YMin) / (YRes - 1)

!	Perform Calculations
	Guess(1) = Theta_Sgnl_Guess
	Guess(2) = Theta_Idlr_Guess
	DO XIndex = 1, XRes, 1
		WRITE(*,"(I3,'% Done')") FLOOR(100.0d0*XIndex/XRes)
		Lambda_Sgnl = XMin + (XIndex - 1)*Lambda_Step 
		Lambda_Idlr=1/(1/Lambda_Pump-1/Lambda_Sgnl)
		Phi_Sgnl = (YMax + YMin) / 2d0
		
		CALL FIND_OPT_PMF_2D(Guess,PMF_Value,Info)

!		Initial guesses for theta (signal,idler) optimum in next iteration
		IF (Use_Previous .AND. (PMF_Value > 0.5)) THEN
       		Guess(1) = Theta_Sgnl
			Guess(2) = Theta_Idlr
		ELSE
			Guess(1) = Theta_Sgnl_Guess
			Guess(2) = Theta_Idlr_Guess
		ENDIF
			
		DO YIndex = 1, YRes, 1
			Phi_Sgnl = YMin + (YIndex - 1)*Phi_Step
			PMF_Value = PHASE_MATCH_FUNC(Success)

			PlotData(XIndex,YIndex,1) = PMF_Value
			XAxis(XIndex) = Lambda_Sgnl
			YAxis(YIndex) = Phi_Sgnl * 180.0d0/Pi
		ENDDO       

		IF ( .NOT. WAVELENGTHS_VALID() ) THEN
			WRITE(*,"('WAVELENGTH VIOLATION - Pump:',F8.4,'; Signal:',F8.4,'; Idler:',F8.4)") Lambda_Pump, Lambda_Sgnl, Lambda_Idlr
			Wavelength_Violation = .TRUE.
		ENDIF
	ENDDO

END SUBROUTINE PLOT_PMF_LAMBDA_PHI

!***************************************************************************
!
!             D Effective Function (Lambda, Phi)
!
!***************************************************************************

SUBROUTINE INFO_DEFF_LAMBDA_PHI
	USE PM_Data

	PlotName = 'D Effective Function (Lambda, Phi)'
	XName = 'Lambda Signal (µm)'
	YName = 'Phi Signal (deg)'
	PlotSize(1) = 50
	PlotSize(2) = 50
	PlotSize(3) = 1
	Plot_Type = PLOT_3D
	XAxisIsAngle = .FALSE.
	YAxisIsAngle = .TRUE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Pump_Width_Pos) = 1
	Params_Used(Crystal_Length_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Theta_Pump_Pos) = 1
	Params_Used(Phi_Pump_Pos) = 1
	Params_Used(TypePM_Pos) = 1
	Params_Used(Use_Previous_Pos) = 1
END SUBROUTINE INFO_DEFF_LAMBDA_PHI

SUBROUTINE PLOT_DEFF_LAMBDA_PHI(XAxis, YAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)
	REAL*8, INTENT(INOUT) :: YAxis(YRes)

!	Local variables
	INTEGER XIndex, YIndex, Info
	REAL*8 PMF_Value, PHASE_MATCH_FUNC, Lambda_Step, Phi_Step, Guess(2), DEff
	REAL*8 Ny, Nz, Argument, GET_D_EFF
	LOGICAL WAVELENGTHS_VALID, Success

!	Set x-axis range
	IF (AutoRange) THEN
		XMin = 1/(1/Lambda_Pump - 1/Abs_Max) 
		XMax = Abs_Max
!		It doesn't like the wavelength right on the edge, so move in a bit
		XMin = XMin+(XMax-XMin)/XRes
		YMin = -0.5d0 * Pi/180d0
		YMax = 0.5d0 * Pi/180d0
		YRes = 50
	END IF
	Lambda_Step = (XMax - XMin) / (XRes - 1)
	Phi_Step = (YMax - YMin) / (YRes - 1)

!	Perform Calculations
	Guess(1) = Theta_Sgnl_Guess
	Guess(2) = Theta_Idlr_Guess
	DO XIndex = 1, XRes, 1
		WRITE(*,"(I3,'% Done')") FLOOR(100.0d0*XIndex/XRes)
		Lambda_Sgnl = XMin + (XIndex - 1)*Lambda_Step 
		Lambda_Idlr=1/(1/Lambda_Pump-1/Lambda_Sgnl)
		IF ( .NOT. WAVELENGTHS_VALID() ) THEN
			WRITE(*,"('WAVELENGTH VIOLATION - Pump:',F8.4,'; Signal:',F8.4,'; Idler:',F8.4)") Lambda_Pump, Lambda_Sgnl, Lambda_Idlr
			Wavelength_Violation = .TRUE.
		ENDIF

		Phi_Sgnl = (YMax + YMin) / 2d0
		Phi_Idlr = Phi_Sgnl + Pi
		
!		CALL FIND_OPT_PMF_2D(Guess,PMF_Value,Info)
		CALL INDEX(Crystal, Lambda_Sgnl, Temperature, N_Sgnl, Ny, Nz)
		CALL INDEX(Crystal, Lambda_Idlr, Temperature, N_Idlr, Ny, Nz)

		Argument = 0.5d0*Lambda_Pump / (N_Pump*N_Sgnl*Lambda_Sgnl) *	&
					( N_Sgnl**2 + (N_Pump*Lambda_Sgnl/Lambda_Pump)**2	&
					 - (N_Idlr*Lambda_Sgnl/Lambda_Idlr)**2 )			
		IF (DABS(Argument) .LE. 1) THEN
			Theta_Sgnl = DACOS(Argument)
		ELSE
			Theta_Sgnl = 0
		ENDIF
							
		Argument = 0.5d0*Lambda_Pump / (N_Pump*N_Idlr*Lambda_Idlr) *	&
					( N_Idlr**2 + (N_Pump*Lambda_Idlr/Lambda_Pump)**2	&
					 - (N_Sgnl*Lambda_Idlr/Lambda_Sgnl)**2 )			
		IF (DABS(Argument) .LE. 1) THEN
			Theta_Idlr = DACOS(Argument)
		ELSE
			Theta_Idlr = 0
		ENDIF

!		Initial guesses for theta (signal,idler) optimum in next iteration
!		IF (Use_Previous .AND. (PMF_Value > 0.5)) THEN
!       		Guess(1) = Theta_Sgnl
!			Guess(2) = Theta_Idlr
!		ELSE
!			Guess(1) = Theta_Sgnl_Guess
!			Guess(2) = Theta_Idlr_Guess
!		ENDIF
			

		DO YIndex = 1, YRes, 1
			Phi_Sgnl = YMin + (YIndex - 1)*Phi_Step
			Phi_Idlr = Phi_Sgnl + Pi

			!PMF_Value = PHASE_MATCH_FUNC(Success)
			CALL Calculate_DEff(DEff)

			PlotData(XIndex,YIndex,1) = DEFF
			XAxis(XIndex) = Lambda_Sgnl
			YAxis(YIndex) = Phi_Sgnl * 180.0d0/Pi
		ENDDO       
	ENDDO
END SUBROUTINE PLOT_DEFF_LAMBDA_PHI

!***************************************************************************
!
!        Angle Spreads Around Optimum Theta Signal (Lambda)
!
!***************************************************************************

SUBROUTINE INFO_PMFSPREAD_LAMBDA
	USE PM_DATA
	IMPLICIT NONE	

!	Begin Subroutine
	PlotName = 'Angle Spreads Around Optimum Theta Signal (Lambda)'
	XName = 'Lambda Signal (µm)'
	YName = 'Theta Signal (deg);Theta Internal Spread FWHM;Theta External Spread FWHM;Internal Gaussian Spread FWHM;Total Internal Spread FWHM'
	PlotSize(1) = 50
	PlotSize(2) = 5
	PlotSize(3) = 1 
	XAxisIsAngle = .FALSE.
	YAxisIsAngle = .FALSE.
	Plot_Type = PLOT_2D
	Params_Used(Pump_Width_Pos) = 1
	Params_Used(Crystal_Length_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Theta_Pump_Pos) = 1
	Params_Used(Phi_Pump_Pos) = 1
	Params_Used(Phi_Sgnl_Pos) = 1
	Params_Used(TypePM_Pos) = 1
	Params_Used(Target_PMF_Pos) = 1
END SUBROUTINE INFO_PMFSPREAD_LAMBDA

SUBROUTINE PLOT_PMFSPREAD_LAMBDA(XAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)

!	Local variables
	INTEGER Info, XIndex
	REAL*8 PMF_Value, Guess(2), Lambda_Step, Argument
	REAL*8 Theta_Sgnl_P, Theta_Sgnl_N, Theta_Sgnl_Ext, Theta_Sgnl_P_Ext, Theta_Sgnl_N_Ext
	REAL*8 Theta_Int_FWHM, Theta_Ext_FWHM, Phi_FWHM, Gaussian_FWHM
	REAL*8 THETA_AT_PMFVALUE, PHI_AT_PMFVALUE
	LOGICAL WAVELENGTHS_VALID

!	Set x-axis range
	IF (AutoRange) THEN
		XMin = 1/(1/Lambda_Pump - 1/Abs_Max) 
		XMax = Abs_Max
!		It doesn't like the wavelength right on the edge, so move in a bit
		XMin = XMin+(XMax-XMin)/XRes
	END IF
	Lambda_Step = (XMax-XMin)/(XRes-1)

!	Perform Calculations
	Guess(1) = Theta_Sgnl_Guess
	Guess(2) = Theta_Idlr_Guess
   
	DO XIndex = 1, XRes, 1
		WRITE(*,"(I3,'% Done')") FLOOR(100.0d0*XIndex/XRes)

		Lambda_Sgnl = XMin + (XIndex - 1)*Lambda_Step 
		Lambda_Idlr = 1/(1/Lambda_Pump - 1/Lambda_Sgnl)

		IF ( .NOT. WAVELENGTHS_VALID() ) THEN
			WRITE(*,"('WAVELENGTH VIOLATION - Pump:',F8.4,'; Signal:',F8.4,'; Idler:',F8.4)") Lambda_Pump, Lambda_Sgnl, Lambda_Idlr
			Wavelength_Violation = .TRUE.
		ENDIF

!		Calculate theta_signal optimum and theta_idler optimum

		CALL FIND_OPT_PMF_2D(Guess,PMF_Value,Info)

!		Give preference to keeping the answer the same sign as the guess
		IF ( (DSIGN(1d0,Theta_Sgnl) .NE. DSIGN(1d0,Guess(1))) .AND. (DSIGN(1d0,Theta_Idlr) .NE. DSIGN(1d0,Guess(2))) ) THEN
			Guess(1) = DSIGN(Theta_Sgnl,Guess(1))
			Guess(2) = DSIGN(Theta_Idlr,Guess(2))
			CALL FIND_OPT_PMF_2D(Guess,PMF_Value,Info)
		ENDIF

!		Initial guesses for theta (signal,idler) optimum in next iteration
		IF ((Use_Previous .NE. 0) .AND. (PMF_Value .GT. 0.05d0)) THEN
       		Guess(1) = Theta_Sgnl
			Guess(2) = Theta_Idlr
		ELSE
			Guess(1) = Theta_Sgnl_Guess
			Guess(2) = Theta_Idlr_Guess
		ENDIF

!		Calculate the spread in Phi_Sgnl
!		Phi_FWHM = DABS(PHI_AT_PMFVALUE(1d0, 0.5d0)-PHI_AT_PMFVALUE(-1d0, 0.5d0))
		Phi_FWHM = 0

!		Calculate the positive and negative theta_signal spread
		Theta_Sgnl_P = THETA_AT_PMFVALUE(1d0, TargetPMFVal)
		Theta_Sgnl_N = THETA_AT_PMFVALUE(-1d0, TargetPMFVal)
		Theta_Int_FWHM = DABS(Theta_Sgnl_P-Theta_Sgnl_N)

!		Convert Internal angles to External angles

		Argument = DSIN( Theta_Sgnl )*N_Sgnl
		IF (DABS(Argument) .LE. 1.0d0) THEN
			Theta_Sgnl_Ext = DASIN( Argument ) 
		ELSE
			IF (Theta_Sgnl .GT. 0d0) THEN
				Theta_Sgnl_Ext = Pi-Theta_Sgnl
			ELSE
				Theta_Sgnl_Ext = -Pi+Theta_Sgnl
			ENDIF						
		ENDIF

		Argument = DSIN(Theta_Sgnl_P)*N_Sgnl
		IF (DABS(Argument) .LE. 1d0) THEN
			Theta_Sgnl_P_Ext = DASIN( Argument )
		ELSE
			Theta_Sgnl_P_Ext = 0d0
		ENDIF

		Argument = DSIN(Theta_Sgnl_N)*N_Sgnl
		IF (DABS(Argument) .LE. 1d0) THEN
			Theta_Sgnl_N_Ext = DASIN( Argument )
		ELSE
			Theta_Sgnl_N_Ext = 0d0
		ENDIF

		IF ((Theta_Sgnl_N_Ext .NE. 0d0) .AND. (Theta_Sgnl_P_Ext .NE. 0d0)) THEN
			Theta_Ext_FWHM = DABS(Theta_Sgnl_P_Ext-Theta_Sgnl_N_Ext)
		ELSE
			Theta_Ext_FWHM = 0d0
		ENDIF

		Gaussian_FWHM = 2*DATAN(Lambda_Sgnl*DLOG(1d0/TargetPMFVal)/Pi/Pump_Width/N_Sgnl/DSQRT(2d0))
!		Write it all out

		XAxis(XIndex) = Lambda_Sgnl
		PlotData(XIndex,1,1) = Theta_Sgnl*180d0/Pi
		PlotData(XIndex,2,1) = Theta_Int_FWHM*180d0/Pi
		PlotData(XIndex,3,1) = Theta_Ext_FWHM*180d0/Pi
		PlotData(XIndex,4,1) = Gaussian_FWHM*180d0/Pi
		PlotData(XIndex,5,1) = DSQRT(Gaussian_FWHM**2+Theta_Int_FWHM**2)*180d0/Pi
	ENDDO
	
END SUBROUTINE PLOT_PMFSPREAD_LAMBDA

!***************************************************************************
!
!              Optimum Theta vs Phi 
!
!***************************************************************************

SUBROUTINE INFO_OPTPMF_THETA_PHI
	USE PM_Data

	PlotName = 'Optimum Thetas (Phi_Sgnl)'
	XName = 'Phi Signal (Degrees)'
	YName = 'Theta Signal;Theta Idler;Theta Signal External;Theta Idler External'
	PlotSize(1) = 722
	PlotSize(2) = 4
	PlotSize(3) = 1
	Plot_Type = PLOT_2D
	XAxisIsAngle = .TRUE.
	YAxisIsAngle = .TRUE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Pump_Width_Pos) = 1
	Params_Used(Crystal_Length_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Theta_Pump_Pos) = 1
	Params_Used(Phi_Pump_Pos) = 1
	Params_Used(TypePM_Pos) = 1
	Params_Used(Lambda_Sgnl_Pos) = 1
END SUBROUTINE INFO_OPTPMF_THETA_PHI

SUBROUTINE PLOT_OPTPMF_THETA_PHI(XAxis, PlotData)
	USE CRYSTAL_DATA; USE PM_DATA
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)

!	Local variables
	INTEGER Info, XIndex, PlotIndex, MaxOnFirst
	REAL*8 Guess(2), Difference, MaxDiff, BestTheta_Sgnl, BestTheta_Idlr
	REAL*8 PMF_Value, Phi_Step, Step
	REAL*8 Argument, Theta_Sgnl_Ext, Theta_Idlr_Ext
	LOGICAL WAVELENGTHS_VALID, MissedPoint

!	Perform Calculations

	IF (AutoRange) THEN
		XMin = 0d0
		XMax = 2d0*Pi
	ENDIF

	Phi_Step = (XMax - XMin) / (FLOOR(XRes/2d0) - 1)
	PlotIndex = 0
	MissedPoint = .FALSE.
	Guess(1) = 0
	Guess(2) = 0

	DO XIndex = 1, FLOOR(XRes/2d0), 1
		WRITE(*,"(I3,'% Done')") FLOOR(100.0d0*XIndex/XRes)
		Phi_Sgnl = XMin + (XIndex - 1) * Phi_Step
		Phi_Idlr = Phi_Sgnl + Pi

		CALL FIND_OPT_PMF_2D(Guess,PMF_Value,Info)

		IF (Use_Previous) THEN
			Guess(1) = Theta_Sgnl
			Guess(2) = Theta_Idlr
		ENDIF

		IF ((PMF_Value .GT. 0.8d0) .AND. (Theta_Sgnl .GT. 0)) THEN
			Argument = N_Sgnl * DSIN(Theta_Sgnl)
			IF (DABS(Argument) .LE. 1.0d0) THEN
				Theta_Sgnl_Ext = DASIN( Argument ) 
			ELSE
				Theta_Sgnl = 0d0
			ENDIF

			Argument = N_Idlr * DSIN(Theta_Idlr)
			IF (DABS(Argument) .LE. 1.0d0) THEN
				Theta_Idlr_Ext = DASIN( Argument ) 
			ELSE
				Theta_Idlr_Ext = 0d0
			ENDIF
	
			PlotIndex = PlotIndex + 1
			XAxis(PlotIndex) = Phi_Sgnl * 180d0 / Pi
			PlotData(PlotIndex,1,1) = Theta_Sgnl * 180d0 / Pi
			PlotData(PlotIndex,2,1) = Theta_Idlr * 180d0 / Pi
			PlotData(PlotIndex,3,1) = Theta_Sgnl_Ext * 180d0 / Pi
			PlotData(PlotIndex,4,1) = Theta_Idlr_Ext * 180d0 / Pi
		ELSE
			MissedPoint = .TRUE.
		ENDIF
	ENDDO

	IF (MissedPoint) THEN
		MaxOnFirst = PlotIndex
		DO XIndex = MaxOnFirst, 1, -1
			WRITE(*,"(I3,'% Done')") FLOOR( 50d0+50d0*(MaxOnFirst-XIndex)/MaxOnFirst)
			Phi_Sgnl = XAxis(XIndex) * Pi/180
			Phi_Idlr = Phi_Sgnl + Pi

			Difference = 0d0
			Step = 0.5d0
			MaxDiff = 0d0
			BestTheta_Sgnl = PlotData(XIndex,1,1) * Pi / 180d0
			BestTheta_Idlr = PlotData(XIndex,2,1) * Pi / 180d0
			Guess(1) = 0
			Guess(2) = 0
			DO WHILE ((MaxDiff .LT. 0.02) .AND. (Guess(1) .LT. (PlotData(XIndex,1,1)+25)*Pi/180d0))
				CALL FIND_OPT_PMF_2D(Guess,PMF_Value,Info)
				Guess(1) = Step * Pi / 180
				Guess(2) = Step * Pi / 180
				Step = Step + 1d0
				Difference = PlotData(XIndex,1,1) * Pi / 180d0 - Theta_Sgnl
				IF (DABS(Difference) .GT. MaxDiff) THEN
					MaxDiff = DABS(Difference)
					BestTheta_Sgnl = Theta_Sgnl
					BestTheta_Idlr = Theta_Idlr
				ENDIF
			END DO
			Theta_Sgnl = BestTheta_Sgnl
			Theta_Idlr = BestTheta_Idlr

			IF ((PMF_Value .GT. 0.8d0) .AND. (Theta_Sgnl .GT. 0)) THEN
				Argument = N_Sgnl * DSIN(Theta_Sgnl)
				IF (DABS(Argument) .LE. 1.0d0) THEN
					Theta_Sgnl_Ext = DASIN( Argument ) 
				ELSE
					Theta_Sgnl = 0d0
				ENDIF

				Argument = N_Idlr * DSIN(Theta_Idlr)
				IF (DABS(Argument) .LE. 1.0d0) THEN
					Theta_Idlr_Ext = DASIN( Argument ) 
				ELSE
					Theta_Idlr_Ext = 0d0
				ENDIF
	
				PlotIndex = PlotIndex + 1
				XAxis(PlotIndex) = Phi_Sgnl * 180d0 / Pi
				PlotData(PlotIndex,1,1) = Theta_Sgnl * 180d0 / Pi
				PlotData(PlotIndex,2,1) = Theta_Idlr * 180d0 / Pi
				PlotData(PlotIndex,3,1) = Theta_Sgnl_Ext * 180d0 / Pi
				PlotData(PlotIndex,4,1) = Theta_Idlr_Ext * 180d0 / Pi
			ENDIF
		ENDDO
	ENDIF

	IF (PlotIndex .GT. 0) Then
		DO XIndex = PlotIndex+1, XRes, 1
			XAxis(XIndex) = XAxis(PlotIndex)
			PlotData(XIndex,1,1) = PlotData(PlotIndex,1,1)
			PlotData(XIndex,2,1) = PlotData(PlotIndex,2,1)
			PlotData(XIndex,3,1) = PlotData(PlotIndex,3,1)
			PlotData(XIndex,4,1) = PlotData(PlotIndex,4,1)
		ENDDO
	ELSE
		WRITE (*,*) 'Unable to find any phasematching for these parameters'
		StatusOK = .FALSE.
	ENDIF
					
	IF ( .NOT. WAVELENGTHS_VALID() ) THEN
		WRITE(*,"('WAVELENGTH VIOLATION - Pump:',F8.4,'; Signal:',F8.4,'; Idler:',F8.4)") Lambda_Pump, Lambda_Sgnl, Lambda_Idlr
		Wavelength_Violation = .TRUE.
	ENDIF
	
END SUBROUTINE PLOT_OPTPMF_THETA_PHI

!***************************************************************************
!
!              Optimum Theta vs Phi 
!
!***************************************************************************

SUBROUTINE INFO_OPTPMF_X_Y
	USE PM_DATA
	IMPLICIT NONE
	
	CALL INFO_OPTPMF_THETA_PHI
	PlotName = 'Optimum Thetas (X,Y)'
	XName = 'X Opening Angle'
	Plot_Type = PLOT_POLAR
END SUBROUTINE INFO_OPTPMF_X_Y

SUBROUTINE PLOT_OPTPMF_X_Y(XAxis, PlotData)
	USE CRYSTAL_DATA; USE PM_DATA
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)

!	Local Variables
	INTEGER I

	CALL PLOT_OPTPMF_THETA_PHI(XAxis, PlotData)
	IF (StatusOK) THEN
		DO I = 1, XRes, 1
			PlotData(I,2,1) = - PlotData(I,2,1)
			PlotData(I,4,1) = - PlotData(I,4,1)
		ENDDO
	ENDIF
END SUBROUTINE PLOT_OPTPMF_X_Y

!***************************************************************************
!
!              Phase Matching function (X,Y)
!
!***************************************************************************

SUBROUTINE INFO_PMF_X_Y
	USE PM_Data

	PlotName = 'Phase Matching function (X,Y)'
	IF (Angle_External) THEN
		XName = 'X Opening Angle, External'
		YName = 'Y Opening Angle, External'
	ELSE
		XName = 'X Opening Angle, Internal'
		YName = 'Y Opening Angle, Internal'
	ENDIF
	PlotSize(1) = 50
	PlotSize(2) = 50
	PlotSize(3) = 1
	Plot_Type = PLOT_3D
	XAxisIsAngle = .TRUE.
	YAxisIsAngle = .TRUE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Pump_Width_Pos) = 1
	Params_Used(Crystal_Length_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Theta_Pump_Pos) = 1
	Params_Used(Phi_Pump_Pos) = 1
	Params_Used(Lambda_Sgnl_Pos) = 1
	Params_Used(TypePM_Pos) = 1
	Params_Used(Angle_External_Pos) = 1
END SUBROUTINE INFO_PMF_X_Y

SUBROUTINE PLOT_PMF_X_Y(XAxis, YAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)
	REAL*8, INTENT(INOUT) :: YAxis(YRes)

!	Local variables
	INTEGER Info, XIndex, YIndex
	REAL*8 Guess, PMF_Value, X, X_Step, Y, Y_Step, Theta_Sgnl_Ext
	LOGICAL WAVELENGTHS_VALID, Success

!	Set axis range
	IF (AutoRange) THEN
		XMin = -5d0 * Pi / 180
		XMax = 5d0 * Pi / 180
		YMin = -5d0 * Pi / 180
		YMax = 5d0 * Pi / 180
	END IF
	X_Step = (XMax - XMin) / (XRes - 1)
	Y_Step = (YMax - YMin) / (YRes - 1)

!	Perform calculations
	Guess = Theta_Idlr_Guess
	
	DO XIndex = 1, XRes, 1
		WRITE(*,"(I3,'% Done')") FLOOR(100.0d0*XIndex/XRes)
		X = DTAN(XMin + (XIndex - 1)*X_Step) 
	
		DO YIndex = 1, YRes, 1
			Y = DTAN(YMin + (YIndex - 1)*Y_Step)
			IF (Angle_External) THEN
				Theta_Sgnl_Ext = DATAN(DSQRT(X**2 + Y**2))
				CALL GetThetaInternal(Theta_Sgnl_Ext,Theta_Sgnl)
			ELSE
				Theta_Sgnl = DATAN(DSQRT(X**2 + Y**2))
			ENDIF
				
			Phi_Sgnl = DATAN(Y/X)
!			Make Phi_Sgnl between 0 and 2*Pi
			IF (X .LT. 0d0) THEN
				Phi_Sgnl = Phi_Sgnl + Pi
			ENDIF
			IF (Phi_Sgnl .LT. 0d0) THEN
				Phi_Sgnl = Phi_Sgnl + 2d0*Pi
			ENDIF
			Phi_Idlr = Phi_Sgnl + Pi

			CALL FIND_OPT_PMF_1D(Guess, PMF_Value, Info)
					
			PlotData(XIndex,YIndex,1) = PlotData(XIndex,YIndex,1)+PMF_Value
			XAxis(XIndex) = X * 180d0 / Pi
			YAxis(YIndex) = Y * 180d0 / Pi
		ENDDO
	ENDDO

	IF ( .NOT. WAVELENGTHS_VALID() ) THEN
		WRITE(*,"('WAVELENGTH VIOLATION - Pump:',F8.4,'; Signal:',F8.4,'; Idler:',F8.4)") Lambda_Pump, Lambda_Sgnl, Lambda_Idlr
		Wavelength_Violation = .TRUE.
	ENDIF
	
END SUBROUTINE PLOT_PMF_X_Y

!***************************************************************************
!
!              Phase Matching function (Theta signal, Theta idler)
!
!***************************************************************************

SUBROUTINE INFO_PMF_THETASIG_THETAIDL
	USE PM_Data

	PlotName = 'Phase Matching function (Theta signal, Theta idler)'
	XName = 'Theta Signal (deg)'
	YName = 'Theta Idler (deg)'
	PlotSize(1) = 100
	PlotSize(2) = 100
	PlotSize(3) = 1
	Plot_Type = PLOT_3D
	XAxisIsAngle = .TRUE.
	YAxisIsAngle = .TRUE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Pump_Width_Pos) = 1
	Params_Used(Crystal_Length_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Lambda_Sgnl_Pos) = 1
	Params_Used(TypePM_Pos) = 1
	Params_Used(Theta_Pump_Pos) = 1
	Params_Used(Phi_Pump_Pos) = 1
	Params_Used(Phi_Sgnl_Pos) = 1
END SUBROUTINE INFO_PMF_THETASIG_THETAIDL

SUBROUTINE PLOT_PMF_THETASIG_THETAIDL(XAxis, YAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)
	REAL*8, INTENT(INOUT) :: YAxis(YRes)

!	Local variables
	INTEGER XIndex, YIndex
	REAL*8 Theta_Sgnl_Step, Theta_Idlr_Step, PHASE_MATCH_FUNC
	LOGICAL WAVELENGTHS_VALID, Success
			        
!	Set x-axis range
	IF (AutoRange) THEN
		XMin = -10d0 * Pi / 180d0
		XMax = 10d0 * Pi / 180d0
		YMin = -10.0d0 * Pi / 180d0
		YMax = 10.0d0 * Pi / 180d0
	END IF
	Theta_Sgnl_Step = (XMax - XMin) / (XRes - 1)
	Theta_Idlr_Step = (YMax - YMin) / (YRes - 1)

!	Perform Calculations
	DO XIndex = 1, XRes, 1
		Theta_Sgnl = XMin + (XIndex - 1)*Theta_Sgnl_Step 
		
		DO YIndex = 1, YRes, 1
			Theta_Idlr = YMin + (YIndex - 1)*Theta_Idlr_Step
		
			PlotData(XIndex,YIndex,1) = PHASE_MATCH_FUNC(Success)
			XAxis(XIndex) = Theta_Sgnl*180/Pi
			YAxis(YIndex) = Theta_Idlr*180/Pi

		END DO
	END DO
	IF ( .NOT. WAVELENGTHS_VALID() ) THEN
		WRITE(*,"('WAVELENGTH VIOLATION - Pump:',F8.4,'; Signal:',F8.4,'; Idler:',F8.4)") Lambda_Pump, Lambda_Sgnl, Lambda_Idlr
		Wavelength_Violation = .TRUE.
	ENDIF

END SUBROUTINE PLOT_PMF_THETASIG_THETAIDL

!***************************************************************************
!
!             Phase Matching function, Signal Direction Fixed (Theta Idler, Phi Idler)
!
!***************************************************************************

SUBROUTINE INFO_PMF_THETAIDL_PHIIDL
	USE PM_Data
	
	PlotName = 'Phase Matching function, Signal Direction Fixed (Theta Idler, Phi Idler)'
	XName = 'Phi Idler (deg)'
	YName = 'Theta Idler (deg)'
	PlotSize(1) = 100
	PlotSize(2) = 100
	PlotSize(3) = 1
	Plot_Type = PLOT_3D
	XAxisIsAngle = .TRUE.
	YAxisIsAngle = .TRUE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Pump_Width_Pos) = 1
	Params_Used(Crystal_Length_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Lambda_Sgnl_Pos) = 1
	Params_Used(TypePM_Pos) = 1
	Params_Used(Theta_Pump_Pos) = 1
	Params_Used(Phi_Pump_Pos) = 1
	Params_Used(Theta_Sgnl_Pos) = 1
	Params_Used(Phi_Sgnl_Pos) = 1
END SUBROUTINE INFO_PMF_THETAIDL_PHIIDL

SUBROUTINE PLOT_PMF_THETAIDL_PHIIDL(XAxis, YAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)
	REAL*8, INTENT(INOUT) :: YAxis(YRes)

!	Local variables
	INTEGER Info, XIndex, YIndex
	REAL*8 Phi_Idlr_Step, Theta_Idlr_Step, PHASE_MATCH_FUNC, PMF_Value, Guess(2)
	REAL*8 Delta_Phi, Delta_Theta, PHI_AT_PMFVALUE, THETA_IDLR_AT_PMFVALUE
	LOGICAL WAVELENGTHS_VALID, Success
	        
!	CALL FIND_OPT_PMF_1D(Theta_Idlr_Guess, PMF_Value, Info)
	IF (AutoRange) THEN
		Guess(1) = 0
		Guess(2) = 0
		CALL FIND_OPT_PMF_2D(Guess,PMF_Value,Info)
		Delta_Phi = DABS(3d0*(PHI_AT_PMFVALUE(1d0,PMF_Value*0.5d0) - Phi_Sgnl))
		IF (Delta_Phi .GT. Pi) THEN
			Delta_Phi = Pi
		ENDIF
		XMin = (Phi_Idlr-Delta_Phi)
		XMax = (Phi_Idlr+Delta_Phi)
		Delta_Theta = DABS(3d0*(THETA_IDLR_AT_PMFVALUE(1d0,PMF_Value*0.5d0) - Theta_Idlr))
		IF (Delta_Phi .GT. Pi/2d0) THEN
			Delta_Phi = Pi/2d0
		ENDIF
		YMin = (Theta_Idlr-Delta_Theta)
		YMax = (Theta_Idlr+Delta_Theta)
	END IF
	Phi_Idlr_Step = (XMax - XMin) / (XRes - 1)
	Theta_Idlr_Step = (YMax - YMin) / (YRes - 1)

!	Perform Calculations
	DO XIndex = 1, XRes, 1
		Phi_Idlr = XMin + (XIndex - 1)*Phi_Idlr_Step 
		
		DO YIndex = 1, YRes, 1
			Theta_Idlr = YMin + (YIndex - 1)*Theta_Idlr_Step

			PlotData(XIndex,YIndex,1) = PHASE_MATCH_FUNC(Success)
			XAxis(XIndex) = Phi_Idlr*180/Pi
			YAxis(YIndex) = Theta_Idlr*180/Pi
		END DO
	END DO
	IF ( .NOT. WAVELENGTHS_VALID() ) THEN
		WRITE(*,"('WAVELENGTH VIOLATION - Pump:',F8.4,'; Signal:',F8.4,'; Idler:',F8.4)") Lambda_Pump, Lambda_Sgnl, Lambda_Idlr
		Wavelength_Violation = .TRUE.
	ENDIF

END SUBROUTINE PLOT_PMF_THETAIDL_PHIIDL


!***************************************************************************
!
!              Optimum Theta Signal at Chosen Phi (Theta_Pump)
!
!***************************************************************************

SUBROUTINE INFO_THETAOPT_THETAPUMP
	USE PM_Data

	PlotName = 'Optimum Theta Signal at Chosen Phi (Theta_Pump)'
	XName = 'Theta Pump (deg)'
	YName = 'Theta Signal Internal (deg);Theta Signal External (deg);Theta of Corellated Photon (deg);Theta External of Correlated Photon(deg);D Effective'
	PlotSize(1) = 50
	PlotSize(2) = 5
	PlotSize(3) = 1
	Plot_Type = PLOT_2D
	XAxisIsAngle = .TRUE.
	YAxisIsAngle = .FALSE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Pump_Width_Pos) = 1
	Params_Used(Crystal_Length_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Lambda_Sgnl_Pos) = 1
	Params_Used(Phi_Pump_Pos) = 1
	Params_Used(Phi_Sgnl_Pos) = 1
	Params_Used(TypePM_Pos) = 1
	Params_Used(Use_Previous_Pos) = 1
END SUBROUTINE INFO_THETAOPT_THETAPUMP

SUBROUTINE PLOT_THETAOPT_THETAPUMP(XAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)

!	Local variables
	INTEGER Info, XIndex
	REAL*8 PMF_Value, Guess(2), NF,NS, Theta_Step, GET_D_EFF, DEff
	REAL*8 Nx, Ny, Nz, Theta_Sgnl_Ext, Theta_Idlr_Ext, Argument
	LOGICAL WAVELENGTHS_VALID
	        
!	Set range
	IF (AutoRange) THEN
		XMin = 0
		XMax = Pi/4d0
	END IF

	Theta_Step = (XMax - XMin) / (XRes - 1)

	Theta_Sgnl_Guess = 0
	Theta_Idlr_Guess = 0

!	Perform calculations
	Guess(1) = Theta_Sgnl_Guess
	Guess(2) = Theta_Idlr_Guess 
	DO XIndex = 1, XRes, 1
		WRITE(*,"(I3,'% Done')") FLOOR(100.0d0*XIndex/XRes)
		Theta_Pump = XMin + (XIndex - 1) * Theta_Step
		CALL Init_NPump


!		Calculate theta_signal optimum and theta_idler optimum
!		IF ( TypePM .EQ. 1) THEN
!
!			CALL INDEX(Crystal, Lambda_Sgnl, Temperature, N_Sgnl, Ny, Nz)
!			CALL INDEX(Crystal, Lambda_Idlr, Temperature, N_Idlr, Ny, Nz)
!
!			Argument = 0.5d0*Lambda_Pump / (N_Pump*N_Sgnl*Lambda_Sgnl) *	&
!						( N_Sgnl**2 + (N_Pump*Lambda_Sgnl/Lambda_Pump)**2	&
!						 - (N_Idlr*Lambda_Sgnl/Lambda_Idlr)**2 )			
!			IF (DABS(Argument) .LE. 1) THEN
!				Theta_Sgnl = DACOS(Argument)
!			ELSE
!				Theta_Sgnl = 0
!			ENDIF
!					
!			Argument = 0.5d0*Lambda_Pump / (N_Pump*N_Idlr*Lambda_Idlr) *	&
!						( N_Idlr**2 + (N_Pump*Lambda_Idlr/Lambda_Pump)**2	&
!						 - (N_Sgnl*Lambda_Idlr/Lambda_Sgnl)**2 )			
!			IF (DABS(Argument) .LE. 1) THEN
!				Theta_Idlr = DACOS(Argument)
!			ELSE
!				Theta_Idlr = 0
!			ENDIF
!		ELSE
		CALL FIND_OPT_PMF_2D(Guess,PMF_Value,Info)

!		Give preference to keeping the answer the same sign as the guess
		IF ( DSIGN(1d0,Theta_Sgnl) .NE. DSIGN(1d0,Theta_Sgnl_Guess) ) THEN
			Guess(1) = DSIGN(Theta_Sgnl,Guess(1))
			Guess(2) = DSIGN(Theta_Idlr,Guess(2))
			CALL FIND_OPT_PMF_2D(Guess,PMF_Value,Info)
		ENDIF

!		Initial guesses for theta (signal,idler) optimum in next iteration
		IF (Use_Previous .AND. (PMF_Value .GT. 0.05d0)) THEN
			Guess(1) = Theta_Sgnl
			Guess(2) = Theta_Idlr
		ELSE
			Guess(1) = Theta_Sgnl_Guess
			Guess(2) = Theta_Idlr_Guess
		ENDIF

		Argument = DSIN( Theta_Sgnl )*N_Sgnl
		IF (DABS(Argument) .LE. 1.0d0) THEN
			Theta_Sgnl_Ext = DASIN( Argument ) 
		ELSE
			IF (Theta_Sgnl .GT. 0d0) THEN
				Theta_Sgnl_Ext = Pi-Theta_Sgnl
			ELSE
				Theta_Sgnl_Ext = -Pi+Theta_Sgnl
			ENDIF						
		ENDIF
		Argument = DSIN( Theta_Idlr )*N_Idlr
		IF (DABS(Argument) .LE. 1.0d0) THEN
			Theta_Idlr_Ext = DASIN( Argument ) 
		ELSE
			IF (Theta_Idlr .GT. 0d0) THEN
				Theta_Idlr_Ext = Pi-Theta_Idlr
			ELSE
				Theta_Idlr_Ext = -Pi+Theta_Idlr
			ENDIF						
		ENDIF

		CALL Calculate_DEff(DEff)
		XAxis(XIndex) = Theta_Pump * 180d0 / Pi
		PlotData(XIndex,1,1) = Theta_Sgnl * 180/Pi
		PlotData(XIndex,2,1) = Theta_Sgnl_Ext * 180/Pi
		PlotData(XIndex,3,1) = Theta_Idlr *180/Pi
		PlotData(XIndex,4,1) = Theta_Idlr_Ext *180d0/Pi
		PlotData(XIndex,5,1) = DEff
	ENDDO
	IF ( .NOT. WAVELENGTHS_VALID() ) THEN
		WRITE(*,"('WAVELENGTH VIOLATION - Pump:',F8.4,'; Signal:',F8.4,'; Idler:',F8.4)") Lambda_Pump, Lambda_Sgnl, Lambda_Idlr
		Wavelength_Violation = .TRUE.
	ENDIF

END SUBROUTINE PLOT_THETAOPT_THETAPUMP

!***************************************************************************
!
!              Phase Match Function (Theta_Pump, Theta_Sgnl)
!
!***************************************************************************

SUBROUTINE INFO_PMF_THETAPUMP_THETASGNL
	USE PM_Data

	PlotName = 'Phase Match Function (Theta_Pump, Theta_Sgnl)'
	XName = 'Theta Pump (deg)'
	YName = 'Theta Signal Internal (deg);Theta Signal External (deg);Theta of Corellated Photon (deg);Theta External of Correlated Photon(deg);D Effective'
	PlotSize(1) = 50
	PlotSize(2) = 50
	PlotSize(3) = 1
	Plot_Type = PLOT_3D
	XAxisIsAngle = .TRUE.
	YAxisIsAngle = .TRUE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Pump_Width_Pos) = 1
	Params_Used(Crystal_Length_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Lambda_Sgnl_Pos) = 1
	Params_Used(Phi_Pump_Pos) = 1
	Params_Used(Phi_Sgnl_Pos) = 1
	Params_Used(TypePM_Pos) = 1
	Params_Used(Use_Previous_Pos) = 1
END SUBROUTINE INFO_PMF_THETAPUMP_THETASGNL

SUBROUTINE PLOT_PMF_THETAPUMP_THETASGNL(XAxis, YAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)
	REAL*8, INTENT(INOUT) :: YAxis(YRes)

!	Local variables
	INTEGER Info, XIndex, YIndex
	REAL*8 PMF_Value, Guess, NF,NS, Theta_P_Step, Theta_S_Step, GET_D_EFF, DEff
	REAL*8 Nx, Ny, Nz, Theta_Sgnl_Ext, Theta_Idlr_Ext, Argument
	LOGICAL WAVELENGTHS_VALID
	        
!	Set range
	IF (AutoRange) THEN
		XMin = 0
		XMax = Pi/4d0
		YMin = 0
		YMax = Pi/4d0
	END IF

	Theta_P_Step = (XMax - XMin) / (XRes - 1)
	Theta_S_Step = (YMax - YMin) / (YRes - 1)

	Theta_Sgnl_Guess = 0
	Theta_Idlr_Guess = 0

!	Perform calculations
	Guess = Theta_Sgnl_Guess
	DO XIndex = 1, XRes, 1
		WRITE(*,"(I3,'% Done')") FLOOR(100.0d0*XIndex/XRes)
		Theta_Pump = XMin + (XIndex - 1) * Theta_P_Step
		CALL Init_NPump

		DO YIndex = 1, YRes, 1
			Theta_Sgnl = YMin + (YIndex - 1) * Theta_S_Step

			CALL FIND_OPT_PMF_1D(Guess,PMF_Value,Info)

!			Initial guesses for next iteration
			IF (Use_Previous .AND. (PMF_Value .GT. 0.05d0)) THEN
				Guess = Theta_Idlr
			ELSE
				Guess = Theta_Idlr_Guess
			ENDIF

			CALL Calculate_DEff(DEff)
			XAxis(XIndex) = Theta_Pump * 180d0 / Pi
			YAxis(YIndex) = Theta_Sgnl * 180d0 / Pi
			PlotData(XIndex,YIndex,1) = PMF_Value
		ENDDO
	ENDDO
	IF ( .NOT. WAVELENGTHS_VALID() ) THEN
		WRITE(*,"('WAVELENGTH VIOLATION - Pump:',F8.4,'; Signal:',F8.4,'; Idler:',F8.4)") Lambda_Pump, Lambda_Sgnl, Lambda_Idlr
		Wavelength_Violation = .TRUE.
	ENDIF

END SUBROUTINE PLOT_PMF_THETAPUMP_THETASGNL

!***************************************************************************
!
!        Angle Spreads Around Optimum Theta Signal (Lambda)
!
!***************************************************************************

SUBROUTINE INFO_PMFSPREAD_THETA_PUMP
	USE PM_DATA
	IMPLICIT NONE	

!	Begin Subroutine
	PlotName = 'Angle Spreads Around Optimum Theta Signal (Theta_Pump)'
	XName = 'Theta Pump (deg)'
	YName = 'Theta Signal (deg);Theta Internal Spread FWHM;Theta External Spread FWHM;Internal Gaussian Spread FWHM;Total Internal Spread FWHM,D Effective at Optimum Theta'
	PlotSize(1) = 50
	PlotSize(2) = 6
	PlotSize(3) = 1 
	XAxisIsAngle = .TRUE.
	YAxisIsAngle = .FALSE.
	Plot_Type = PLOT_2D
	Params_Used(Pump_Width_Pos) = 1
	Params_Used(Crystal_Length_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Lambda_Sgnl_Pos) = 1
	Params_Used(Phi_Pump_Pos) = 1
	Params_Used(Phi_Sgnl_Pos) = 1
	Params_Used(TypePM_Pos) = 1
	Params_Used(Target_PMF_Pos) = 1
END SUBROUTINE INFO_PMFSPREAD_THETA_PUMP

SUBROUTINE PLOT_PMFSPREAD_THETA_PUMP(XAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)

!	Local variables
	INTEGER Info, XIndex
	REAL*8 PMF_Value, Guess(2), Theta_Step, Argument
	REAL*8 Theta_Sgnl_P, Theta_Sgnl_N, Theta_Sgnl_Ext, Theta_Sgnl_P_Ext, Theta_Sgnl_N_Ext
	REAL*8 Theta_Int_FWHM, Theta_Ext_FWHM, Phi_FWHM, Gaussian_FWHM
	REAL*8 THETA_AT_PMFVALUE, PHI_AT_PMFVALUE, DEff
	LOGICAL WAVELENGTHS_VALID

!	Set x-axis range
	IF (AutoRange) THEN
		XMin = 0 
		XMax = Pi / 4
	END IF
	Theta_Step = (XMax-XMin)/(XRes-1)

!	Perform Calculations
	Guess(1) = Theta_Sgnl_Guess
	Guess(2) = Theta_Idlr_Guess
   
	DO XIndex = 1, XRes, 1
		WRITE(*,"(I3,'% Done')") FLOOR(100.0d0*XIndex/XRes)
		Theta_Pump = XMin + (XIndex - 1)*Theta_Step 
		CALL Init_NPump

		IF ( .NOT. WAVELENGTHS_VALID() ) THEN
			WRITE(*,"('WAVELENGTH VIOLATION - Pump:',F8.4,'; Signal:',F8.4,'; Idler:',F8.4)") Lambda_Pump, Lambda_Sgnl, Lambda_Idlr
			Wavelength_Violation = .TRUE.
		ENDIF

!		Calculate theta_signal optimum and theta_idler optimum

		CALL FIND_OPT_PMF_2D(Guess,PMF_Value,Info)

!		Give preference to keeping the answer the same sign as the guess
		IF ( (DSIGN(1d0,Theta_Sgnl) .NE. DSIGN(1d0,Guess(1))) .AND. (DSIGN(1d0,Theta_Idlr) .NE. DSIGN(1d0,Guess(2))) ) THEN
			Guess(1) = DSIGN(Theta_Sgnl,Guess(1))
			Guess(2) = DSIGN(Theta_Idlr,Guess(2))
			CALL FIND_OPT_PMF_2D(Guess,PMF_Value,Info)
		ENDIF

!		Initial guesses for theta (signal,idler) optimum in next iteration
		IF ((Use_Previous .NE. 0) .AND. (PMF_Value .GT. 0.05d0)) THEN
       		Guess(1) = Theta_Sgnl
			Guess(2) = Theta_Idlr
		ELSE
			Guess(1) = Theta_Sgnl_Guess
			Guess(2) = Theta_Idlr_Guess
		ENDIF

		CALL Calculate_DEff(DEff)

!		Calculate the spread in Phi_Sgnl
!		Phi_FWHM = DABS(PHI_AT_PMFVALUE(1d0, 0.5d0)-PHI_AT_PMFVALUE(-1d0, 0.5d0))
		Phi_FWHM = 0

!		Calculate the positive and negative theta_signal spread
		Theta_Sgnl_P = THETA_AT_PMFVALUE(1d0, TargetPMFVal)
		Theta_Sgnl_N = THETA_AT_PMFVALUE(-1d0, TargetPMFVal)
		Theta_Int_FWHM = DABS(Theta_Sgnl_P-Theta_Sgnl_N)

!		Convert Internal angles to External angles

		Argument = DSIN( Theta_Sgnl )*N_Sgnl
		IF (DABS(Argument) .LE. 1.0d0) THEN
			Theta_Sgnl_Ext = DASIN( Argument ) 
		ELSE
			IF (Theta_Sgnl .GT. 0d0) THEN
				Theta_Sgnl_Ext = Pi-Theta_Sgnl
			ELSE
				Theta_Sgnl_Ext = -Pi+Theta_Sgnl
			ENDIF						
		ENDIF

		Argument = DSIN(Theta_Sgnl_P)*N_Sgnl
		IF (DABS(Argument) .LE. 1d0) THEN
			Theta_Sgnl_P_Ext = DASIN( Argument )
		ELSE
			Theta_Sgnl_P_Ext = 0d0
		ENDIF

		Argument = DSIN(Theta_Sgnl_N)*N_Sgnl
		IF (DABS(Argument) .LE. 1d0) THEN
			Theta_Sgnl_N_Ext = DASIN( Argument )
		ELSE
			Theta_Sgnl_N_Ext = 0d0
		ENDIF

		IF ((Theta_Sgnl_N_Ext .NE. 0d0) .AND. (Theta_Sgnl_P_Ext .NE. 0d0)) THEN
			Theta_Ext_FWHM = DABS(Theta_Sgnl_P_Ext-Theta_Sgnl_N_Ext)
		ELSE
			Theta_Ext_FWHM = 0d0
		ENDIF

		Gaussian_FWHM = 2*DATAN(Lambda_Sgnl*DLOG(1d0/TargetPMFVal)/Pi/Pump_Width/N_Sgnl/DSQRT(2d0))
!		Write it all out

		XAxis(XIndex) = Theta_Pump*180d0/Pi
		PlotData(XIndex,1,1) = Theta_Sgnl*180d0/Pi
		PlotData(XIndex,2,1) = Theta_Int_FWHM*180d0/Pi
		PlotData(XIndex,3,1) = Theta_Ext_FWHM*180d0/Pi
		PlotData(XIndex,4,1) = Gaussian_FWHM*180d0/Pi
		PlotData(XIndex,5,1) = DSQRT(Gaussian_FWHM**2+Theta_Int_FWHM**2)*180d0/Pi
		PlotData(XIndex,6,1) = DEff
	ENDDO
	
END SUBROUTINE PLOT_PMFSPREAD_THETA_PUMP

!***************************************************************************
!
!              ANGLE_TEST
!
!***************************************************************************

SUBROUTINE INFO_ANGLE_TEST
	USE PM_Data

	PlotName = 'Angle Test'
	XName = 'Vertical Tilt (deg)'
	YName = 'Theta Tilt;Theta Tilt Ext;Phi Tilt;Theta Pump;Phi Pump;X;Y;Z'
	PlotSize(1) = 50
	PlotSize(2) = 8
	PlotSize(3) = 1
	Plot_Type = PLOT_2D
	XAxisIsAngle = .TRUE.
	YAxisIsAngle = .FALSE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Theta_Cut_Pos) = 1
	Params_Used(Phi_Cut_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Cube_Tilt_H_Pos) = 1
	Params_Used(Cube_Tilt_V_Pos) = 1
END SUBROUTINE INFO_ANGLE_TEST

SUBROUTINE PLOT_ANGLE_TEST(XAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)

!	Local variables
	INTEGER XIndex
	REAL*8 Tilt_Step
	        
!	Set range
	IF (AutoRange) THEN
		XMin = 0
		XMax = Pi/4d0
	END IF

	Tilt_Step = (XMax - XMin) / (XRes - 1)

	DO XIndex = 1, XRes, 1
		WRITE(*,"(I3,'% Done')") FLOOR(100.0d0*XIndex/XRes)
		Cube_Tilt_V = XMin + (XIndex - 1) * Tilt_Step

		CALL SET_INTERNAL_ANGLES

		XAxis(XIndex) = Cube_Tilt_V * 180d0 / Pi
		PlotData(XIndex,1,1) = Theta_Tilt * 180d0 / Pi
		PlotData(XIndex,2,1) = Theta_Tilt_Ext * 180d0 / Pi
		PlotData(XIndex,3,1) = Phi_Tilt * 180d0 / Pi
		PlotData(XIndex,4,1) = Theta_Pump * 180d0 / Pi
		PlotData(XIndex,5,1) = Phi_Pump * 180d0 / Pi
		PlotData(XIndex,6,1) = DSIN(Theta_Tilt) * DCOS(Phi_Tilt)
		PlotData(XIndex,7,1) = DSIN(Theta_Tilt) * DSIN(Phi_Tilt)
		PlotData(XIndex,8,1) = DCOS(Theta_Tilt) 


	ENDDO

END SUBROUTINE PLOT_ANGLE_TEST

!***************************************************************************
!
!              Phase Match Function (Theta_Pump, Theta_Sgnl)
!
!***************************************************************************

SUBROUTINE INFO_REFLECTANCE_PHI
	USE PM_Data

	PlotName = 'Reflectance (Phi Signal)'
	XName = 'Phi Signal (deg)'
	YName = 'Signal Reflectance;Idler Reflectance'
	PlotSize(1) = 50
	PlotSize(2) = 2
	PlotSize(3) = 1
	Plot_Type = PLOT_2D
	XAxisIsAngle = .TRUE.
	YAxisIsAngle = .FALSE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Pump_Width_Pos) = 1
	Params_Used(Crystal_Length_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Lambda_Sgnl_Pos) = 1
	Params_Used(Theta_Pump_Pos) = 1
	Params_Used(Phi_Pump_Pos) = 1
	Params_Used(TypePM_Pos) = 1
	Params_Used(Use_Previous_Pos) = 1
	Params_Used(Cube_Tilt_H_Pos) = 1
	Params_Used(Cube_Tilt_V_Pos) = 1
	Params_Used(Theta_Cut_Pos) = 1
	Params_Used(Phi_Cut_Pos) = 1
END SUBROUTINE INFO_REFLECTANCE_PHI

SUBROUTINE PLOT_REFLECTANCE_PHI(XAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)

!	Local variables
	INTEGER Info, XIndex
	REAL*8 Phi_Step, Guess(2)
	REAL*8 Ref_Sgnl, Ref_Idlr, PMF_Value
	LOGICAL WAVELENGTHS_VALID

!	Set range
	IF (AutoRange) THEN
		XMin = 0
		XMax = 2d0*Pi
	END IF

	Phi_Step = (XMax - XMin) / (XRes - 1)
	Theta_Idlr_Guess = 0

	CALL SET_INTERNAL_ANGLES
	CALL INIT_NPump

	Guess(1) = 0d0
	Guess(2) = 0d0
	CALL FIND_OPT_PMF_2D(Guess,PMF_Value,Info)

!	Perform calculations
	DO XIndex = 1, XRes, 1
		WRITE(*,"(I3,'% Done')") FLOOR(100.0d0*XIndex/XRes)
		Phi_Sgnl = XMin + (XIndex - 1) * Phi_Step
		Phi_Idlr = Phi_Sgnl + Pi

		CALL GET_REFLECTANCE(Lambda_Sgnl,Theta_Sgnl,Phi_Sgnl,0,Ref_Sgnl)
		CALL GET_REFLECTANCE(Lambda_Idlr,Theta_Idlr,Phi_Idlr,0,Ref_Idlr)

		XAxis(XIndex) = Phi_Sgnl * 180d0 / Pi
		PlotData(XIndex,1,1) = Ref_Sgnl
		PlotData(XIndex,2,1) = Ref_Idlr

	ENDDO
	IF ( .NOT. WAVELENGTHS_VALID() ) THEN
		WRITE(*,"('WAVELENGTH VIOLATION - Pump:',F8.4,'; Signal:',F8.4,'; Idler:',F8.4)") Lambda_Pump, Lambda_Sgnl, Lambda_Idlr
		Wavelength_Violation = .TRUE.
	ENDIF

END SUBROUTINE PLOT_REFLECTANCE_PHI

!***************************************************************************
!
!              Phase Match Function (Theta_Pump, Theta_Sgnl)
!
!***************************************************************************

SUBROUTINE INFO_REFLECTANCE_LAMBDA_PHI
	USE PM_Data

	PlotName = 'Signal Reflectance (Phi Signal)'
	XName = 'Lambda Signal (um)'
	YName = 'Phi Signal'
	PlotSize(1) = 50
	PlotSize(2) = 50
	PlotSize(3) = 1
	Plot_Type = PLOT_3D
	XAxisIsAngle = .FALSE.
	YAxisIsAngle = .TRUE.
	Params_Used(Crystal_Pos) = 1
	Params_Used(Pump_Width_Pos) = 1
	Params_Used(Crystal_Length_Pos) = 1
	Params_Used(Temperature_Pos) = 1
	Params_Used(Lambda_Pump_Pos) = 1
	Params_Used(Theta_Pump_Pos) = 1
	Params_Used(Phi_Pump_Pos) = 1
	Params_Used(TypePM_Pos) = 1
	Params_Used(Use_Previous_Pos) = 1
	Params_Used(Cube_Tilt_H_Pos) = 1
	Params_Used(Cube_Tilt_V_Pos) = 1
	Params_Used(Theta_Cut_Pos) = 1
	Params_Used(Phi_Cut_Pos) = 1
END SUBROUTINE INFO_REFLECTANCE_LAMBDA_PHI

SUBROUTINE PLOT_REFLECTANCE_LAMBDA_PHI(XAxis, YAxis, PlotData)
	USE PM_Data; USE Crystal_Data
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(INOUT) :: PlotData(XRes,YRes,ZRes)
	REAL*8, INTENT(INOUT) :: XAxis(XRes)
	REAL*8, INTENT(INOUT) :: YAxis(YRes)

!	Local variables
	INTEGER Info, XIndex, YIndex
	REAL*8 Phi_Step, Lambda_Step, Guess(2)
	REAL*8 Ref_Sgnl, Ref_Idlr, PMF_Value
	LOGICAL WAVELENGTHS_VALID

!	Set range
	IF (AutoRange) THEN
		YMin = 0
		YMax = 2d0*Pi
		XMin = 0.699d0
		XMax = 0.705d0
	END IF

	Phi_Step = (YMax - YMin) / (YRes - 1)
	Lambda_Step = (XMax - XMin) / (XRes - 1)
	Theta_Idlr_Guess = 0

	CALL SET_INTERNAL_ANGLES
	CALL INIT_NPump

	Guess(1) = 0d0
	Guess(2) = 0d0

!	Perform calculations
	DO XIndex = 1, XRes, 1
		WRITE(*,"(I3,'% Done')") FLOOR(100.0d0*XIndex/XRes)
		Lambda_Sgnl = XMin + (XIndex - 1) * Lambda_Step

		CALL FIND_OPT_PMF_2D(Guess,PMF_Value,Info)
		Guess(1) = Theta_Sgnl
		Guess(2) = Theta_Idlr

		DO YIndex = 1, YRes, 1
			Phi_Sgnl = YMin + (YIndex - 1) * Phi_Step
			Phi_Idlr = Phi_Sgnl + Pi


			CALL GET_REFLECTANCE(Lambda_Sgnl,Theta_Sgnl,Phi_Sgnl,0,Ref_Sgnl)
			CALL GET_REFLECTANCE(Lambda_Idlr,Theta_Idlr,Phi_Idlr,0,Ref_Idlr)

			XAxis(XIndex) = Lambda_Sgnl
			YAxis(YIndex) = Phi_Sgnl * 180d0 / Pi
			PlotData(XIndex,YIndex,1) = Ref_Sgnl

		ENDDO
	ENDDO
	IF ( .NOT. WAVELENGTHS_VALID() ) THEN
		WRITE(*,"('WAVELENGTH VIOLATION - Pump:',F8.4,'; Signal:',F8.4,'; Idler:',F8.4)") Lambda_Pump, Lambda_Sgnl, Lambda_Idlr
		Wavelength_Violation = .TRUE.
	ENDIF

END SUBROUTINE PLOT_REFLECTANCE_LAMBDA_PHI