!**************************************************************************************
!     
!  TRANSFORM_CUBE_TO_CRYSTAL(CubeVector, CrystalVector)
!                                                                                     
!***************************************************************************************

SUBROUTINE TRANSFORM_CUBE_TO_CRYSTAL(CubeVector, CrystalVector)
	USE PM_Data
	IMPLICIT NONE

!	Argument
	
	REAL*8, INTENT(IN)  :: CubeVector(3)
	REAL*8, INTENT(OUT) :: CrystalVector(3)

!	Local Variables

	REAL*8 Cos_Theta_Cut, Sin_Theta_Cut, Cos_Phi_Cut, Sin_Phi_Cut

!	Perform Transformation

	Cos_Theta_Cut = DCOS(Theta_Cut)
	Sin_Theta_Cut = DSIN(Theta_Cut)
	Cos_Phi_Cut   = DCOS(Phi_Cut)
	Sin_Phi_Cut   = DSIN(Phi_Cut)

	CrystalVector(1) =   Cos_Theta_Cut * Cos_Phi_Cut * CubeVector(1) &
				       -                 Sin_Phi_Cut * CubeVector(2) &
				       + Sin_Theta_Cut * Cos_Phi_Cut * CubeVector(3)

	CrystalVector(2) =   Cos_Theta_Cut * Sin_Phi_Cut * CubeVector(1) &
				       +                 Cos_Phi_Cut * CubeVector(2) &
				       + Sin_Theta_Cut * Sin_Phi_Cut * CubeVector(3)

	CrystalVector(3) = - Sin_Theta_Cut * CubeVector(1) &
		               + Cos_Theta_Cut * CubeVector(3)

END SUBROUTINE TRANSFORM_CUBE_TO_CRYSTAL

!**************************************************************************************
!     
!  TRANSFORM_CUBE_TO_CRYSTAL(CubeVector, CrystalVector)
!                                                                                     
!***************************************************************************************

SUBROUTINE TRANSFORM_CUBE_TO_TABLE(CubeVector, TableVector)
	USE PM_Data
	IMPLICIT NONE

!	Argument
	
	REAL*8, INTENT(IN)  :: CubeVector(3)
	REAL*8, INTENT(OUT) :: TableVector(3)

!	Local Variables

	REAL*8 Cos_H_Tilt, Sin_H_Tilt, Cos_V_Tilt, Sin_V_Tilt

!	Perform Transformation

	Cos_H_Tilt = DCOS(Cube_Tilt_H)
	Sin_H_Tilt = DSIN(Cube_Tilt_H)
	Cos_V_Tilt = DCOS(Cube_Tilt_V)
	Sin_V_Tilt = DSIN(Cube_Tilt_V)

	TableVector(1) =   Cos_H_Tilt              * CubeVector(1) &
	                 + Sin_H_Tilt * Sin_V_Tilt * CubeVector(2) &
	                 + Sin_H_Tilt * Cos_V_Tilt * CubeVector(2)


END SUBROUTINE TRANSFORM_CUBE_TO_TABLE

!**************************************************************************************
!     
!  TRANSFORM_PUMP_TO_CUBE(PumpVector, CubeVector)
!                                                                                     
!***************************************************************************************

SUBROUTINE TRANSFORM_PUMP_TO_CUBE(PumpVector, CubeVector)
	USE PM_Data
	IMPLICIT NONE

!	Argument
	
	REAL*8, INTENT(IN)  :: PumpVector(3)
	REAL*8, INTENT(OUT) :: CubeVector(3)

!	Local Variables

	REAL*8 Cos_Theta_Tilt, Sin_Theta_Tilt, Cos_Phi_Tilt, Sin_Phi_Tilt

!	Perform Transformation

	Cos_Theta_Tilt = DCOS(Theta_Tilt)
	Sin_Theta_Tilt = DSIN(Theta_Tilt)
	Cos_Phi_Tilt   = DCOS(Phi_Tilt)
	Sin_Phi_Tilt   = DSIN(Phi_Tilt)

	CubeVector(1) =   Cos_Theta_Tilt * Cos_Phi_Tilt * PumpVector(1) &
				    -                 Sin_Phi_Tilt * PumpVector(2) &
	                + Sin_Theta_Tilt * Cos_Phi_Tilt * PumpVector(3)

	CubeVector(2) =   Cos_Theta_Tilt * Sin_Phi_Tilt * PumpVector(1) &
	                +                 Cos_Phi_Tilt * PumpVector(2) &
	                + Sin_Theta_Tilt * Sin_Phi_Tilt * PumpVector(3)

	CubeVector(3) = - Sin_Theta_Tilt * PumpVector(1) &
	                + Cos_Theta_Tilt * PumpVector(3)

END SUBROUTINE TRANSFORM_Pump_TO_Cube

!**************************************************************************************
!     
!  TRANSFORM_PUMP_TO_CRYSTAL(PumpVector, CubeVector)
!                                                                                     
!***************************************************************************************

SUBROUTINE TRANSFORM_PUMP_TO_CRYSTAL(PumpVector, CrystalVector)
	USE PM_Data
	IMPLICIT NONE

!	Argument
	
	REAL*8, INTENT(IN)  :: PumpVector(3)
	REAL*8, INTENT(OUT) :: CrystalVector(3)

!	Local Variables

	REAL*8 Cos_Theta_Pump, Sin_Theta_Pump, Cos_Phi_Pump, Sin_Phi_Pump

!	Perform Transformation

	Cos_Theta_Pump = DCOS(Theta_Pump)
	Sin_Theta_Pump = DSIN(Theta_Pump)
	Cos_Phi_Pump   = DCOS(Phi_Pump)
	Sin_Phi_Pump   = DSIN(Phi_Pump)

	CrystalVector(1) =   Cos_Theta_Pump * Cos_Phi_Pump * PumpVector(1) &
				       -                  Sin_Phi_Pump * PumpVector(2) &
	                   + Sin_Theta_Pump * Cos_Phi_Pump * PumpVector(3)

	CrystalVector(2) =   Cos_Theta_Pump * Sin_Phi_Pump * PumpVector(1) &
	                   +                  Cos_Phi_Pump * PumpVector(2) &
	                   + Sin_Theta_Pump * Sin_Phi_Pump * PumpVector(3)

	CrystalVector(3) = - Sin_Theta_Pump * PumpVector(1) &
	                   + Cos_Theta_Pump * PumpVector(3)

END SUBROUTINE TRANSFORM_PUMP_TO_CRYSTAL

!***************************************************************************
!
!      SUBROUTINE GET_INDICES(NF, NS, Crystal, Lambda, Theta, Phi)
!
!       Computes the fast and slow indices for a given direction 
!		(Angles measured in crystal coordinates)
!
!		ARGUMENTS	
!		NF      = OUTPUT: Fast index of refraction
!		NS      = OUTPUT: Slow index of refraction
!		Crystal = INPUT: Crystal number
!		Lambda  = INPUT: Wavelength (microns)
!		Theta   = INPUT: Polar angle with respect to z-axis
!		Phi     = INPUT: Azymuthal angle about z-axis with respect to x-axis
!
!***************************************************************************
SUBROUTINE GET_INDICES(NF, NS, Crystal, Lambda, Temperature, Theta, Phi)
	IMPLICIT NONE

	INTEGER Crystal
	REAL*8 NF, NS, Lambda, Theta, Phi, Temperature
	REAL*8 Sx, Sy, Sz, B, C
	REAL*8 Nx,Ny, Nz, ABS_MIN, ABS_MAX

!	Compute the components of the S unit vector
	Sx = DCOS(Phi)*DSIN(Theta)
	Sy = DSIN(Phi)*DSIN(Theta)
	Sz = DCOS(Theta)

!	Get the indices for the Crystal
	CALL INDEX(Crystal, Lambda, Temperature, Nx, Ny, Nz)

!	Compute the two indices in the S direction
	B = (Sx**2)*(1/Ny**2 + 1/Nz**2) + &
		(Sy**2)*(1/Nx**2 + 1/Nz**2) + &
		(Sz**2)*(1/Nx**2 + 1/Ny**2)
	C = (Sx**2)/(Ny**2 * Nz**2) + &
		(Sy**2)/(Nx**2 * Nz**2) + &
		(Sz**2)/(Nx**2 * Ny**2)

	NF = DSQRT( 2.0d0/(B + DSQRT(B**2 - 4*C)) )
	NS = DSQRT( 2.0d0/(B - DSQRT(B**2 - 4*C)) )
END SUBROUTINE GET_INDICES

!***************************************************************************
!
!      Subroutine: INIT_NPUMP
!
!      Initializes N_Pump to the Fast index
!
!***************************************************************************
SUBROUTINE INIT_NPUMP
	USE PM_DATA
	IMPLICIT NONE

!	Local variables

	REAL*8 NF, NS      

!	Begin the subroutine

	CALL GET_INDICES(NF, NS, Crystal, Lambda_Pump, Temperature, Theta_Pump, Phi_Pump)
	N_Pump = NF

END SUBROUTINE INIT_NPUMP

!***************************************************************************
!
!      CALC_INDICES(Lambda, S, NF, NS)
!
!       Computes the fast and slow indices for a given direction 
!
!		ARGUMENTS	
!		NF      = OUTPUT: Fast index of refraction
!		NS      = OUTPUT: Slow index of refraction
!		Lambda  = INPUT: Wavelength (microns)
!		S       = INPUT: Unit vector in direction of interest(In Crystal Coords)
!
!***************************************************************************
SUBROUTINE CALC_INDICES(Lambda, S, NF, NS)
	USE PM_Data
	IMPLICIT NONE

	REAL*8, INTENT(IN)  :: Lambda, S(3)
	REAL*8, INTENT(OUT) :: NF, NS

!	Local Variables
	REAL*8 Nx, Ny, Nz, B, C

!	Compute the two indices in the S direction

	CALL INDEX(Crystal, Lambda, Temperature, Nx, Ny, Nz)

	B = (S(1)**2)*(1/Ny**2 + 1/Nz**2) + &
		(S(2)**2)*(1/Nx**2 + 1/Nz**2) + &
		(S(3)**2)*(1/Nx**2 + 1/Ny**2)

	C = (S(1)**2)/(Ny**2 * Nz**2) + &
		(S(2)**2)/(Nx**2 * Nz**2) + &
		(S(3)**2)/(Nx**2 * Ny**2)

	NF = DSQRT( 2d0/(B + DSQRT(B**2 - 4d0*C)) )
	NS = DSQRT( 2d0/(B - DSQRT(B**2 - 4d0*C)) )

END SUBROUTINE CALC_INDICES

!**************************************************************************************
!     
! GET_REFLECTANCE(Theta,Phi,Pol): Calculate reflectance for a beam at Theta and Phi in the 
!                                 pump frame, with ordinary (Pol=0) or extraordinary (Pol<>0)
!                                 polarization
!                                                                                     
!***************************************************************************************

SUBROUTINE GET_REFLECTANCE(Lambda,Theta,Phi,Pol,Reflectance)
	USE PM_DATA; USE CRYSTAL_DATA
	IMPLICIT NONE

	REAL*8, INTENT(IN)  :: Lambda,Theta, Phi
	INTEGER, INTENT(IN) :: Pol
	REAL*8, INTENT(OUT) :: Reflectance

!	Local Variables
	REAL*8 S_Pump(3), S_Cube(3), S_Crystal(3)
	REAL*8 Po_Crystal(3), Pe_Crystal(3)
	REAL*8 Ns(3), Np(3), Cube_Normal(3), Norm, Es, Ep
	REAL*8 I_In, I_Out, N_Fast, N_Slow, Index, Arg, rs, rp

!	Get the direction vector in the pump frame

	S_Pump(1) = DSIN(Theta) * DCOS(Phi)
	S_Pump(2) = DSIN(Theta) * DSIN(Phi)
	S_Pump(3) = DCOS(Theta) 
	
!	Transform to the Crystal frame

	CALL TRANSFORM_PUMP_TO_CRYSTAL(S_Pump,S_Crystal)

!	Get the cube z-vector (surface normal) in crystal coords

	Cube_Normal(1) = DSIN(Theta_Cut) * DCOS(Phi_Cut)
	Cube_Normal(2) = DSIN(Theta_Cut) * DSIN(Phi_Cut)
	Cube_Normal(3) = DCOS(Theta_Cut) 

!	Get the S and P polarization directions

	Ns(1) = S_Crystal(2) * Cube_Normal(3) - S_Crystal(3) * Cube_Normal(2)
	Ns(2) = S_Crystal(3) * Cube_Normal(1) - S_Crystal(1) * Cube_Normal(3)
	Ns(3) = S_Crystal(1) * Cube_Normal(2) - S_Crystal(2) * Cube_Normal(1)
	Norm = DSQRT( Ns(1)**2 + Ns(2)**2 + Ns(3)**2 )
	Ns(1) = Ns(1) / Norm
	Ns(2) = Ns(2) / Norm
	Ns(3) = Ns(3) / Norm

	Np(1) = S_Crystal(2) * Ns(3) - S_Crystal(3) * Ns(2)
	Np(2) = S_Crystal(3) * Ns(1) - S_Crystal(1) * Ns(3)
	Np(3) = S_Crystal(1) * Ns(2) - S_Crystal(2) * Ns(1)


!	Get the ordinary and extraordinary polarization direction in the crystal frame

	Po_Crystal(1) =  S_Crystal(2) / DSQRT( S_Crystal(1)**2 + S_Crystal(2)**2 )
	Po_Crystal(2) = -S_Crystal(1) / DSQRT( S_Crystal(1)**2 + S_Crystal(2)**2 )
	Po_Crystal(3) = 0d0

	Pe_Crystal(1) = - S_Crystal(3) * Po_Crystal(2)
	Pe_Crystal(2) =   S_Crystal(3) * Po_Crystal(1)
	Pe_Crystal(3) =   S_Crystal(1) * Po_Crystal(2) - S_Crystal(2) * Po_Crystal(1)

	IF (Pol .EQ. 0) THEN
		Es = DABS( Po_Crystal(1)*Ns(1) + Po_Crystal(2)*Ns(2) + Po_Crystal(3)*Ns(3) )
		Ep = DABS( Po_Crystal(1)*Np(1) + Po_Crystal(2)*Np(2) + Po_Crystal(3)*Np(3) )
	ELSE
		Es = DABS( Pe_Crystal(1)*Ns(1) + Pe_Crystal(2)*Ns(2) + Pe_Crystal(3)*Ns(3) )
		Ep = DABS( Pe_Crystal(1)*Np(1) + Pe_Crystal(2)*Np(2) + Pe_Crystal(3)*Np(3) )
	ENDIF

!	Get the index of refraction
	
	CALL CALC_INDICES(Lambda, S_Crystal, N_Fast, N_Slow)
	IF (CrystalType .EQ. NegativeUniaxial) THEN
		IF (Pol .EQ. 0) THEN
			Index = N_Slow
		ELSE
			Index = N_Fast
		ENDIF
	ELSE
		IF (Pol .NE. 0) THEN
			Index = N_Fast
		ELSE
			Index = N_Slow
		ENDIF
	ENDIF

!	Get the internal angle of incidence

	I_In = DACOS(DABS(  S_Crystal(1)*Cube_Normal(1) + S_Crystal(2)*Cube_Normal(2) &
	                  + S_Crystal(3)*Cube_Normal(3) ))

!	Check for internal reflection and calculate reflectance
	Arg = Index * DSIN(I_In)
	IF (Arg .LE. 1) THEN
		I_Out = DASIN(Index * DSIN(I_In) )
		rs = (Index * DCOS(I_In) - DCOS(I_Out))/(Index * DCOS(I_In) + DCOS(I_Out))
		rp = (Index * DCOS(I_Out) - DCOS(I_In))/(Index * DCOS(I_Out) + DCOS(I_In))
		Reflectance = (rs*Es)**2 + (rp*Ep)**2
	ELSE
		Reflectance = 1d0
	ENDIF

!	
END SUBROUTINE GET_REFLECTANCE

!**************************************************************************************
!     
! SET_INTERNAL_ANGLES : Set the internal angles given a Cube Tilt
!                                                                                     
!***************************************************************************************

SUBROUTINE SET_INTERNAL_ANGLES
	USE PM_DATA; USE CRYSTAL_DATA
	IMPLICIT NONE

!	Local Variables
	REAL*8 Sx, Sy, S_Cube(3), S_Crystal(3)
	REAL*8 Cos_Phi_Tilt, Sin_Phi_Tilt, Err
	INTEGER Info

	EXTERNAL Theta_Pump_Int_Err

!	Get Phi Tilt
	IF (Cube_Tilt_V .EQ. 0d0) THEN
		Phi_Tilt = 0d0
	ELSE
		Sx = -DSIN(Cube_Tilt_H)
		Sy =  DCOS(Cube_Tilt_H) * DSIN(Cube_Tilt_V)
		Phi_Tilt = DACOS(Sx / DSQRT(Sx**2 + Sy**2))
		IF (Sy .LT. 0) Phi_Tilt = 2d0 * Pi - Phi_Tilt
	ENDIF

!	Get the external angle of incidence
	Theta_Tilt_Ext = DASIN(DSQRT(DCOS(Cube_Tilt_H)**2 * DSIN(Cube_Tilt_V)**2 + DSIN(Cube_Tilt_H)**2 ))

!	Get the internal angle of incidence
	CALL MINIMIZE_1D(Theta_Tilt_Ext,Theta_Pump_Int_Err,Theta_Tilt,Err,Info)

!	Get Pump Direction in crystal frame
	S_Cube(1) = DSIN(Theta_Tilt) * DCOS(Phi_Tilt)
	S_Cube(2) = DSIN(Theta_Tilt) * DSIN(Phi_Tilt)
	S_Cube(3) = DCOS(Theta_Tilt) 
	CALL TRANSFORM_CUBE_TO_CRYSTAL(S_Cube,S_Crystal)

	Phi_Pump = DACOS( S_Crystal(1) / DSQRT( S_Crystal(1)**2 + S_Crystal(2)**2 ) )
	IF (S_Crystal(2) .LT. 0) Phi_Pump = 2d0 * Pi - Phi_Pump

	Theta_Pump = DASIN(DSQRT( S_Crystal(1)**2 + S_Crystal(2)**2 ))
	IF (S_Crystal(3) .LT. 0) Theta_Pump = Pi - Theta_Pump
	
END SUBROUTINE SET_INTERNAL_ANGLES

!**************************************************************************************
!     
! Theta_Pump_Int_Err(Theta_Param, Err) : Error in snells law for a given Theta_int
!                                                                                     
!***************************************************************************************

SUBROUTINE Theta_Pump_Int_Err(Theta_Param, Err)
	USE PM_DATA
	IMPLICIT NONE

	REAL*8 Theta_Param, Err, Get_N_Sgnl, S_Cube(3), S_Crystal(3), NF, NS
	LOGICAL Success

!	Get direction in cube frame
	
	Theta_Tilt = Theta_Param
	S_Cube(1) = DSIN(Theta_Tilt) * DCOS(Phi_Tilt)
	S_Cube(2) = DSIN(Theta_Tilt) * DSIN(Phi_Tilt)
	S_Cube(3) = DCOS(Theta_Tilt) 

!	Transform to Crystal frame

	CALL TRANSFORM_CUBE_TO_CRYSTAL(S_Cube,S_Crystal)
	CALL CALC_INDICES(Lambda_Pump, S_Crystal, NF, NS)

!	Calculate error with snells law

	Err = DABS( NF * DSIN(Theta_Tilt) - DSIN(Theta_Tilt_Ext) )

END SUBROUTINE Theta_Pump_Int_Err

!**************************************************************************************
!     
!  GetThetaInternal : Given an external Theta_Sgnl, we want to find to 
!                     corresponding internal Theta_Sgnl.  The following 
!                     three subroutines do this
!                                                                                     
!***************************************************************************************

SUBROUTINE GetThetaInternal(Theta_Sgnl_External, Theta_Sgnl_Internal)
	USE PM_DATA
	IMPLICIT NONE

	REAL*8, INTENT(IN) :: Theta_Sgnl_External
	REAL*8, INTENT(OUT) :: Theta_Sgnl_Internal
	COMMON /FMINSHARE/ ThetaExternal
	REAL*8 ThetaExternal
	REAL*8 F
	INTEGER Info
	EXTERNAL FMIN

	ThetaExternal = Theta_Sgnl_External
	CALL MINIMIZE_1D(Theta_Sgnl,FMIN,Theta_Sgnl,F,Info)
	Theta_Sgnl_Internal = Theta_Sgnl
END SUBROUTINE GetThetaInternal

SUBROUTINE FMIN(Theta, F)
	USE PM_DATA
	IMPLICIT NONE

	REAL*8 Theta, F, Get_N_Sgnl
	LOGICAL Success

	COMMON /FMINSHARE/ ThetaExternal
	REAL*8 ThetaExternal

	Theta_Sgnl = Theta
	F = DABS(Get_N_Sgnl(Success) * DSIN(Theta_Sgnl) - DSIN(ThetaExternal))
END SUBROUTINE FMIN

REAL*8 FUNCTION Get_N_Sgnl(Success)
	USE PM_DATA
	IMPLICIT NONE

!	Local variables
       
	LOGICAL Success
	REAL*8 Nx_Sgnl, Ny_Sgnl, Nz_Sgnl
	REAL*8 Sx_Sgnl, Sy_Sgnl, Sz_Sgnl
	REAL*8 B_Sgnl, C_Sgnl

			
	Success = .True.

!	Get Sx, Sy, and Sz in crystal coordinates for signal

	Sx_Sgnl =   DCOS(Phi_Pump) * DCOS(Theta_Pump) * DSIN(Theta_Sgnl) * DCOS(Phi_Sgnl) &
			  + DCOS(Phi_Pump) * DSIN(Theta_Pump) * DCOS(Theta_Sgnl) &
			  - DSIN(Phi_Pump) * DSIN(Theta_Sgnl) * DSIN(Phi_Sgnl)

	Sy_Sgnl =   DSIN(Phi_Pump) * DCOS(Theta_Pump) * DSIN(Theta_Sgnl) * DCOS(Phi_Sgnl) &
			  + DSIN(Phi_Pump) * DSIN(Theta_Pump) * DCOS(Theta_Sgnl) &
			  + DCOS(Phi_Pump) * DSIN(Theta_Sgnl) * DSIN(Phi_Sgnl)

	Sz_Sgnl =   DCOS(Theta_Pump) * DCOS(Theta_Sgnl) &
			  - DSIN(Theta_Pump) * DSIN(Theta_Sgnl) * DCOS(Phi_Sgnl)

!	Get the crystal indices of refraction for the signal wavelength

	CALL INDEX(Crystal, Lambda_Sgnl, Temperature, Nx_Sgnl,Ny_Sgnl, Nz_Sgnl)

!	Calculate the B and C coefficents for signal

	B_Sgnl = (Sx_Sgnl**2) * (1/Ny_Sgnl**2 + 1/Nz_Sgnl**2) + &
			 (Sy_Sgnl**2) * (1/Nx_Sgnl**2 + 1/Nz_Sgnl**2) +	&
			 (Sz_Sgnl**2) * (1/Nx_Sgnl**2 + 1/Ny_Sgnl**2)

	C_Sgnl = (Sx_Sgnl**2)/(Ny_Sgnl**2 * Nz_Sgnl**2) +	&
			 (Sy_Sgnl**2)/(Nx_Sgnl**2 * Nz_Sgnl**2) +	&
			 (Sz_Sgnl**2)/(Nx_Sgnl**2 * Ny_Sgnl**2)

!	Get the signal index (fast or slow depending on phase matching type)

	IF ( TypePM .EQ. 1 ) THEN
		Get_N_Sgnl = DSQRT( 2.0d0/(B_Sgnl - DSQRT(B_Sgnl**2 - 4*C_Sgnl)) )
	ELSE IF ( TypePM .EQ. 2 ) THEN
		Get_N_Sgnl = DSQRT( 2.0d0/(B_Sgnl + DSQRT(B_Sgnl**2 - 4*C_Sgnl)) )
	ELSE
		Get_N_Sgnl = DSQRT( 2.0d0/(B_Sgnl - DSQRT(B_Sgnl**2 - 4*C_Sgnl)) )
	ENDIF

END FUNCTION Get_N_Sgnl
!**************************************************************************************
!      
!  SUBROUTINE Calculate_DeltaK(N,Theta,F)                                                       
!                                                                                     
!	Calculates delta k for parameters in PM_Data module
!	Returns Success = .False. on an error
!
!***************************************************************************************
SUBROUTINE Calculate_DEff(DEff)
	USE PM_DATA; USE CRYSTAL_DATA
	IMPLICIT NONE

!	Arguments

	REAL*8, INTENT(OUT) :: DEff

!	Local variables
	INTEGER I,J
	REAL*8 Sx_Idlr, Sx_Sgnl, Sy_Idlr, Sy_Sgnl, Sz_Idlr, Sz_Sgnl
	REAL*8 Px_Idlr, Px_Sgnl, Py_Idlr, Py_Sgnl, Pz_Idlr, Pz_Sgnl
	REAL*8 P_Pump(3)
	REAL*8 Norm, F(6)

!	Get Sx, Sy, and Sz in crystal coordinates for signal

	IF ((CrystalType .EQ. PositiveUniaxial) .OR. (CrystalType .EQ. NegativeUniaxial)) THEN

		Sx_Sgnl =   DCOS(Theta_Pump) * DCOS(Phi_Pump) * DSIN(Theta_Sgnl) * DCOS(Phi_Sgnl) &
				  -                    DSIN(Phi_Pump) * DSIN(Theta_Sgnl) * DSIN(Phi_Sgnl) &
				  + DSIN(Theta_Pump) * DCOS(Phi_Pump) * DCOS(Theta_Sgnl)

		Sy_Sgnl =   DCOS(Theta_Pump) * DSIN(Phi_Pump) * DSIN(Theta_Sgnl) * DCOS(Phi_Sgnl) &
				  +                    DCOS(Phi_Pump) * DSIN(Theta_Sgnl) * DSIN(Phi_Sgnl) &
				  + DSIN(Theta_Pump) * DSIN(Phi_Pump) * DCOS(Theta_Sgnl)

		Sz_Sgnl = - DSIN(Theta_Pump)                  * DSIN(Theta_Sgnl) * DCOS(Phi_Sgnl) &
		          + DCOS(Theta_Pump)                  * DCOS(Theta_Sgnl) 

!		This should be fixed to account for the possibility of a positive crystal
		Norm = DSQRT( Sx_Sgnl**2 + Sy_Sgnl**2 )
		IF ((TypePM .EQ. 1) .OR. (TypePM .EQ. 3)) THEN
!			Ordinary
			Px_Sgnl = Sy_Sgnl / Norm
			Py_Sgnl = -Sx_Sgnl / Norm
			Pz_Sgnl = 0
		ELSE
!			Extraordinary
			Px_Sgnl = - Sx_Sgnl * Sz_Sgnl / Norm
			Py_Sgnl = Sy_Sgnl * Sz_Sgnl / Norm
			Pz_Sgnl = (Sx_Sgnl**2 - Sy_Sgnl**2) / Norm
		ENDIF

!		Get Sx, Sy, and Sz in crystal coordinates for idler

		Sx_Idlr =   DCOS(Theta_Pump) * DCOS(Phi_Pump) * DSIN(Theta_Idlr) * DCOS(Phi_Idlr) &
				  -                    DSIN(Phi_Pump) * DSIN(Theta_Idlr) * DSIN(Phi_Idlr) &
				  + DSIN(Theta_Pump) * DCOS(Phi_Pump) * DCOS(Theta_Idlr)

		Sy_Idlr =   DCOS(Theta_Pump) * DSIN(Phi_Pump) * DSIN(Theta_Idlr) * DCOS(Phi_Idlr) &
				  +                    DCOS(Phi_Pump) * DSIN(Theta_Idlr) * DSIN(Phi_Idlr) &
				  + DSIN(Theta_Pump) * DSIN(Phi_Pump) * DCOS(Theta_Idlr)

		Sz_Idlr = - DSIN(Theta_Pump)                  * DSIN(Theta_Idlr) * DCOS(Phi_Idlr) &
		          + DCOS(Theta_Pump)                  * DCOS(Theta_Idlr) 


!		This should be fixed to account for the possibility of a positive crystal
		Norm = DSQRT( Sx_Idlr**2 + Sy_Idlr**2 )
		IF ((TypePM .EQ. 1) .OR. (TypePM .EQ. 2)) THEN
!			Ordinary
			Px_Idlr = Sy_Idlr / Norm
			Py_Idlr = -Sx_Idlr / Norm
			Pz_Idlr = 0
		ELSE
!			Extraordinary
			Px_Idlr = - Sx_Idlr * Sz_Idlr / Norm
			Py_Idlr = Sy_Idlr * Sz_Idlr / Norm
			Pz_Idlr = (Sx_Idlr**2 - Sy_Idlr**2) / Norm
		ENDIF

!		This should be fixed to account for the possibility of a positive crystal
!		P_Pump(1) = DSIN(Phi_Pump)
!		P_Pump(2) = -DCOS(Phi_Pump)
!		P_Pump(3) = 0

		P_Pump(1) = DCOS(Theta_Pump) * DCOS(Phi_Pump)
		P_Pump(2) = DCOS(Theta_Pump) * DSIN(Phi_Pump)
		P_Pump(3) = -DSIN(Theta_Pump)
		
		F(1) = Px_Sgnl * Px_Idlr
		F(2) = Py_Sgnl * Py_Idlr
		F(3) = Pz_Sgnl * Pz_Idlr
		F(4) = Py_Sgnl * Pz_Idlr + Pz_Sgnl * Py_Idlr
		F(5) = Pz_Sgnl * Px_Idlr + Px_Sgnl * Pz_Idlr
		F(6) = Px_Sgnl * Py_Idlr + Py_Sgnl * Px_Idlr

		DEff = 0
		DO I = 1, 3, 1
			DO J = 1, 6, 1
				DEff = Deff + D(I,J) * F(J) * P_Pump(I)
			ENDDO
		ENDDO
		DEff = DABS(DEff)

	ENDIF
END SUBROUTINE Calculate_DEff


!**************************************************************************************
!      
!  SUBROUTINE Calculate_DeltaK(N,Theta,F)                                                       
!                                                                                     
!	Calculates delta k for parameters in PM_Data module
!	Returns Success = .False. on an error
!
!***************************************************************************************
SUBROUTINE Calculate_DeltaK
	USE PM_DATA
	IMPLICIT NONE

!	Local variables
	REAL*8 Nx_Idlr, Nx_Sgnl, Ny_Idlr, Ny_Sgnl, Nz_Idlr, Nz_Sgnl
	REAL*8 Sx_Idlr, Sx_Sgnl, Sy_Idlr, Sy_Sgnl, Sz_Idlr, Sz_Sgnl
	REAL*8 B_Idlr, B_Sgnl, C_Idlr, C_Sgnl

!	Get Sx, Sy, and Sz in crystal coordinates for signal

	Sx_Sgnl =   DCOS(Theta_Pump) * DCOS(Phi_Pump) * DSIN(Theta_Sgnl) * DCOS(Phi_Sgnl) &
			  -                    DSIN(Phi_Pump) * DSIN(Theta_Sgnl) * DSIN(Phi_Sgnl) &
			  + DSIN(Theta_Pump) * DCOS(Phi_Pump) * DCOS(Theta_Sgnl) 

	Sy_Sgnl =   DCOS(Theta_Pump) * DSIN(Phi_Pump) * DSIN(Theta_Sgnl) * DCOS(Phi_Sgnl) &
			  +                    DCOS(Phi_Pump) * DSIN(Theta_Sgnl) * DSIN(Phi_Sgnl) &
			  + DSIN(Theta_Pump) * DSIN(Phi_Pump) * DCOS(Theta_Sgnl) 

	Sz_Sgnl = - DSIN(Theta_Pump)                  * DSIN(Theta_Sgnl) * DCOS(Phi_Sgnl) &
	          + DCOS(Theta_Pump)                  * DCOS(Theta_Sgnl)
			  

!	Get the index of refraction for the signal

	CALL INDEX(Crystal, Lambda_Sgnl, Temperature, Nx_Sgnl,Ny_Sgnl, Nz_Sgnl)

!	Calculate the B and C coefficents for signal

	B_Sgnl = (Sx_Sgnl**2) * (1/Ny_Sgnl**2 + 1/Nz_Sgnl**2) + &
			 (Sy_Sgnl**2) * (1/Nx_Sgnl**2 + 1/Nz_Sgnl**2) +	&
			 (Sz_Sgnl**2) * (1/Nx_Sgnl**2 + 1/Ny_Sgnl**2)

	C_Sgnl = (Sx_Sgnl**2)/(Ny_Sgnl**2 * Nz_Sgnl**2) +	&
			 (Sy_Sgnl**2)/(Nx_Sgnl**2 * Nz_Sgnl**2) +	&
			 (Sz_Sgnl**2)/(Nx_Sgnl**2 * Ny_Sgnl**2)

!	Get Sx, Sy, and Sz in crystal coordinates for idler

	Sx_Idlr =   DCOS(Theta_Pump) * DCOS(Phi_Pump) * DSIN(Theta_Idlr) * DCOS(Phi_Idlr) &
			  -                    DSIN(Phi_Pump) * DSIN(Theta_Idlr) * DSIN(Phi_Idlr) &
			  + DSIN(Theta_Pump) * DCOS(Phi_Pump) * DCOS(Theta_Idlr) 

	Sy_Idlr =   DCOS(Theta_Pump) * DSIN(Phi_Pump) * DSIN(Theta_Idlr) * DCOS(Phi_Idlr) &
			  +                    DCOS(Phi_Pump) * DSIN(Theta_Idlr) * DSIN(Phi_Idlr) &
			  + DSIN(Theta_Pump) * DSIN(Phi_Pump) * DCOS(Theta_Idlr) 

	Sz_Idlr = -	DSIN(Theta_Pump)                  * DSIN(Theta_Idlr) * DCOS(Phi_Idlr) &
	          + DCOS(Theta_Pump)                  * DCOS(Theta_Idlr) 
			  

!	Get the index of refraction for the idler

	CALL INDEX(Crystal, Lambda_Idlr, Temperature, Nx_Idlr, Ny_Idlr, Nz_Idlr)

!	Calculate the B and C coefficents for idler

	B_Idlr = (Sx_Idlr**2)*(1/Ny_Idlr**2 + 1/Nz_Idlr**2) + &
			 (Sy_Idlr**2)*(1/Nx_Idlr**2 + 1/Nz_Idlr**2) + &
			 (Sz_Idlr**2)*(1/Nx_Idlr**2 + 1/Ny_Idlr**2)

	C_Idlr = (Sx_Idlr**2)/(Ny_Idlr**2 * Nz_Idlr**2) + &
			 (Sy_Idlr**2)/(Nx_Idlr**2 * Nz_Idlr**2) + &
			 (Sz_Idlr**2)/(Nx_Idlr**2 * Ny_Idlr**2)

!	Get the signal and idler indices (fast or slow depending on phase matching type)

	IF ( TypePM .EQ. 1 ) THEN
              
		N_Sgnl = DSQRT( 2.0d0/(B_Sgnl - DSQRT(B_Sgnl**2 - 4*C_Sgnl)) )
		N_Idlr = DSQRT( 2.0d0/(B_Idlr - DSQRT(B_Idlr**2 - 4*C_Idlr)) )
              
	ELSE IF ( TypePM .EQ. 2 ) THEN
              
		N_Sgnl = DSQRT( 2.0d0/(B_Sgnl + DSQRT(B_Sgnl**2 - 4*C_Sgnl)) )
		N_Idlr = DSQRT( 2.0d0/(B_Idlr - DSQRT(B_Idlr**2 - 4*C_Idlr)) )

	ELSE IF ( TypePM .EQ. 3 ) THEN

		N_Sgnl = DSQRT( 2.0d0/(B_Sgnl - DSQRT(B_Sgnl**2 - 4*C_Sgnl)) )
		N_Idlr = DSQRT( 2.0d0/(B_Idlr + DSQRT(B_Idlr**2 - 4*C_Idlr)) )

	ENDIF

!	Get Sx, Sy, and Sz in pump (") coordinates so we can comput Delta K

	Sx_Sgnl = DSIN(Theta_Sgnl)*DCOS(Phi_Sgnl)
	Sy_Sgnl = DSIN(Theta_Sgnl)*DSIN(Phi_Sgnl)
	Sz_Sgnl = DCOS(Theta_Sgnl)

	Sx_Idlr = DSIN(Theta_Idlr)*DCOS(Phi_Idlr)
	Sy_Idlr = DSIN(Theta_Idlr)*DSIN(Phi_Idlr)
	Sz_Idlr = DCOS(Theta_Idlr)

!	Calculate the components of Delta K

	DKxyz(1) = 2*Pi *(N_Sgnl*Sx_Sgnl/Lambda_Sgnl + N_Idlr*Sx_Idlr/Lambda_Idlr)
	DKxyz(2) = 2*Pi *(N_Sgnl*Sy_Sgnl/Lambda_Sgnl + N_Idlr*Sy_Idlr/Lambda_Idlr)
	DKxyz(3) = 2*Pi *(N_Sgnl*Sz_Sgnl/Lambda_Sgnl + N_Idlr*Sz_Idlr/Lambda_Idlr	  &
				- N_Pump/Lambda_Pump)

END SUBROUTINE Calculate_DeltaK

!***************************************************************************
!
!      FUNCTION: PHASE_MATCH_FUNC(Success)
!
!		Computes the phase matching function for current values of 
!		the global variables in PM_Data.  Success = .False. on a weird result
!
!***************************************************************************
REAL*8 FUNCTION PHASE_MATCH_FUNC(Success)
	USE PM_DATA
	IMPLICIT NONE

	LOGICAL Success

!	Local variables
	REAL*8 SincArg

!	Begin the function

	CALL Calculate_DeltaK
	SincArg = 0.5d0 * Crystal_Length * DKxyz(3)
	IF (SincArg .EQ. 0.0d0) THEN
		PHASE_MATCH_FUNC = DEXP( (-0.5d0) * Pump_Width**2 * (DKxyz(1)**2 + DKxyz(2)**2) )
	ELSE
		PHASE_MATCH_FUNC = DEXP( (-0.5d0) * Pump_Width**2 * (DKxyz(1)**2 + DKxyz(2)**2) ) * (DSIN(SincArg)/SincArg)**2  
	ENDIF
	IF ((PHASE_MATCH_FUNC .GE. 0d0) .AND. (PHASE_MATCH_FUNC .LE. 1d0)) THEN
		Success = .TRUE.
	ELSE
		Success = .FALSE.
	ENDIF
END FUNCTION PHASE_MATCH_FUNC

!**************************************************************************************
!      
!  SUBROUTINE FIND_OPT_PMF_1D(Guess, PMF_Value, Info)
!  Find the optimum phase match function by letting Theta_Idlr vary
!
!***************************************************************************************

SUBROUTINE FIND_OPT_PMF_1D(Guess, PMF_Value, Info)
	USE PM_DATA
	IMPLICIT NONE

!	Arguments

	REAL*8, INTENT(INOUT) :: Guess
	REAL*8, INTENT(OUT) :: PMF_Value
	INTEGER, INTENT(OUT) :: Info

!	Local variables
	REAL*8 Value, PHASE_MATCH_FUNC
	LOGICAL Success
	EXTERNAL PMFSmooth_1D, PMFSearch_1D

!	Begin Subroutine

!	Get close using the smooth function minimization
	CALL MINIMIZE_1D(Guess,PMFSmooth_1D,Theta_Idlr,Value,Info)

!	If we are not in the range of validity for the smooth function,
!	we have to keep looking with the real PMF

	IF (DABS(Crystal_Length*DKxyz(3)) .LT. Pi) THEN
		Guess = Theta_Idlr
		CALL MINIMIZE_1D(Guess,PMFSearch_1D,Theta_Idlr,Value,Info)
	ENDIF

	PMF_Value = PHASE_MATCH_FUNC(Success)
END SUBROUTINE FIND_OPT_PMF_1D

!**************************************************************************************
!      
!  SUBROUTINE PMFSmooth_1D(Theta,F)
!  Calculates negative Log of Phase Match Function and smooths out sinc ripples for search
!  algorythm
!
!***************************************************************************************
SUBROUTINE PMFSmooth_1D(Theta_Idler_Param,Value)
	USE PM_DATA
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(IN) :: Theta_Idler_Param
	REAL*8, INTENT(OUT) :: Value

!	Local variables
	REAL*8 SincArg
	LOGICAL Success

!	Begin the function

	Theta_Idlr = Theta_Idler_Param
	CALL Calculate_DeltaK
	SincArg = 0.5d0 * Crystal_Length * DKxyz(3)
	IF (SincArg .EQ. 0.0d0) THEN
		Value = 0.5d0 * Pump_Width**2 * (DKxyz(1)**2 + DKxyz(2)**2) 
	ELSEIF (DABS(SincArg) .GT. Pi/2.0d0) THEN
		Value = 0.5d0 * Pump_Width**2 * (DKxyz(1)**2 + DKxyz(2)**2) + DLOG(SincArg**2)
	ELSE
		Value = 0.5d0 * Pump_Width**2 * (DKxyz(1)**2 + DKxyz(2)**2) - DLOG((DSIN(SincArg)/SincArg)**2)
	ENDIF

END SUBROUTINE PMFSmooth_1D

!**************************************************************************************
!      
!  SUBROUTINE PMFSearch_1D(Theta,F)                                                       
!  Calculates negative Log of Phase Match Function with only Theta_Idler as a parameter
!
!***************************************************************************************
SUBROUTINE PMFSearch_1D(Theta_Idler_Param,Value)
	USE PM_DATA
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(IN) :: Theta_Idler_Param
	REAL*8, INTENT(OUT) :: Value

!	Local variables
	REAL*8 SincArg
	LOGICAL Success

!	Begin the function

	Theta_Idlr = Theta_Idler_Param
	CALL Calculate_DeltaK
	SincArg = 0.5d0 * Crystal_Length * DKxyz(3)
	IF (SincArg .EQ. 0.0d0) THEN
		Value = 0.5d0 * Pump_Width**2 * (DKxyz(1)**2 + DKxyz(2)**2) 
	ELSE
		IF (DSIN(SincArg) .EQ. 0d0) THEN
			SincArg = SincArg + 1d-11
		ENDIF
		Value = 0.5d0 * Pump_Width**2 * (DKxyz(1)**2 + DKxyz(2)**2) - DLOG((DSIN(SincArg)/SincArg)**2)
	ENDIF

END SUBROUTINE PMFSearch_1D

!**************************************************************************************
!      
!  SUBROUTINE FIND_OPT_PMF_2D(Guess, PMF_Value, Info)
!  Find the optimum phase match function by letting both
!  Theta_Sgnl and Theta_Idlr vary
!
!***************************************************************************************
SUBROUTINE FIND_OPT_PMF_2D(Guess, PMF_Value, Info)
	USE PM_DATA
	IMPLICIT NONE

!	Arguments

	REAL*8, INTENT(INOUT) :: Guess(2)
	REAL*8, INTENT(OUT) :: PMF_Value
	INTEGER, INTENT(OUT) :: Info

!	Local variables
	REAL*8 F, Theta(2), PHASE_MATCH_FUNC
	LOGICAL Success
	EXTERNAL PMFSmooth_2D, PMFSearch_2D

!	Begin Subroutine

!	Get close using the smooth function minimization
	CALL MINIMIZE_2D(Guess,PMFSmooth_2D,Theta,F,Info)

!	If we are not in the range of validity for the smooth function,
!	we have to keep looking with the real PMF

	IF (DABS(Crystal_Length*DKxyz(3)) .LT. Pi) THEN
		Guess(1) = Theta_Sgnl
		Guess(2) = Theta_Idlr
		CALL MINIMIZE_2D(Guess,PMFSearch_2D,Theta,F,Info)
	ENDIF

	Theta_Sgnl = Theta(1)
	Theta_Idlr = Theta(2)

	PMF_Value = PHASE_MATCH_FUNC(Success)
END SUBROUTINE FIND_OPT_PMF_2D

!**************************************************************************************
!      
!  SUBROUTINE PMFSearch_2D(Theta,F)                                                       
!                                                                                     
!  Calculates negative Log of Phase Match Function with both
!  Theta_Sgnl and Theta_Idler as a parameters
!
!***************************************************************************************
SUBROUTINE PMFSearch_2D(Theta,Value)
	USE PM_DATA
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(IN) :: Theta(2)
	REAL*8, INTENT(OUT) :: Value

!	Local variables
	REAL*8 SincArg
	LOGICAL Success

!	Begin the function

	Theta_Sgnl = Theta(1)
	Theta_Idlr = Theta(2)
	CALL Calculate_DeltaK
	SincArg = 0.5d0 * Crystal_Length * DKxyz(3)
	IF (SincArg .EQ. 0.0d0) THEN
		Value = 0.5d0 * Pump_Width**2 * (DKxyz(1)**2 + DKxyz(2)**2) 
	ELSE
		IF (DSIN(SincArg) .EQ. 0d0) THEN
			SincArg = SincArg + 1d-11
		ENDIF
		Value = 0.5d0 * Pump_Width**2 * (DKxyz(1)**2 + DKxyz(2)**2) - DLOG((DSIN(SincArg)/SincArg)**2)
	ENDIF

END SUBROUTINE PMFSearch_2D


!**************************************************************************************
!      
!  SUBROUTINE PMFSmooth_2D(Theta,F)                                                       
!                                                                                     
!  Calculates negative log of the Phase Match Function and replaces sine with 1 when
!  its argument is greater than Pi/2.  This smooths out sinc ripples and increases 
!  relief for the search algorythm. Inputs both theta_Sgnl and Theta_Idlr
!
!***************************************************************************************
SUBROUTINE PMFSmooth_2D(Theta,Value)
	USE PM_DATA
	IMPLICIT NONE

!	Arguments
	REAL*8, INTENT(IN) :: Theta(2)
	REAL*8, INTENT(OUT) :: Value

!	Local variables
	REAL*8 SincArg
	LOGICAL Success

!	Begin the function
	Theta_Sgnl = Theta(1)
	Theta_Idlr = Theta(2)
	CALL Calculate_DeltaK
	SincArg = 0.5d0 * Crystal_Length * DKxyz(3)
	IF (SincArg .EQ. 0.0d0) THEN
		Value = 0.5d0 * Pump_Width**2 * (DKxyz(1)**2 + DKxyz(2)**2) 
	ELSEIF (DABS(SincArg) .GT. Pi/2.0d0) THEN
		Value = 0.5d0 * Pump_Width**2 * (DKxyz(1)**2 + DKxyz(2)**2) + DLOG(SincArg**2)
	ELSE
		Value = 0.5d0 * Pump_Width**2 * (DKxyz(1)**2 + DKxyz(2)**2) - DLOG((DSIN(SincArg)/SincArg)**2)
	ENDIF

END SUBROUTINE PMFSmooth_2D




!***************************************************************************
!
! THETA_Sgnl_AT_PMFVALUE : Looks for a Theta_Sgnl where the phase matching function 
!                          matches a specified value.  Starts looking in the direction
!                          (+,-) specified by the argument Direction.
!
!***************************************************************************
REAL*8 FUNCTION THETA_AT_PMFVALUE(Direction, TargetPMF)
	USE PM_DATA
	IMPLICIT NONE

!	Arguments

	REAL*8, INTENT(IN) :: Direction
	REAL*8, INTENT(IN) :: TargetPMF

!	Local Variables

	REAL*8 X1, X2, F1, F2, Temp, PHASE_MATCH_FUNC, XLow, XHigh, FOut, FIn, XStart, XStop
	REAL*8 Error, FNew, XNew, Theta_Sgnl_In, Max_Theta_Out, XStep, XTest, XAnswer, Theta_Idlr_In
	INTEGER Info, Iterations, Peaks
	LOGICAL Success, InRange, Done

!	Begin Function
	XAnswer = Theta_Sgnl
	Theta_Sgnl_In = Theta_Sgnl
	Theta_Idlr_In = Theta_Idlr

	X1 = Theta_Sgnl
	X2 = Theta_Sgnl + DSIGN(Pi/2d0,Direction)
!	Give preference to finding a value of the same sign as input Theta_Sgnl
	IF ((DSIGN(1d0,X2) .NE. DSIGN(1d0,X1))) THEN
		X2 = 0d0
		Theta_Sgnl = X2
		CALL FIND_OPT_PMF_1D(Theta_Idlr_Guess,F2,Info)
		F2 = F2 - TargetPMF
		IF (F2 .GT. 0d0) THEN
			X2 = Theta_Sgnl + DSIGN(Pi/2d0,Direction)
		ENDIF
	ENDIF

!	Make sure we have a point where the PMF Value is bigger than target
	Theta_Sgnl = X1
	CALL FIND_OPT_PMF_1D(Theta_Idlr_Guess,F1,Info)
	F1 = F1 - TargetPMF
	IF (F1 .GT. 0d0) THEN

!		Do a Binary Search to Bracket at least one zero crossing
		Done = .FALSE.
		Iterations = 0
		DO WHILE ((.NOT. Done) .AND. (Iterations .LT. 400))
			XNew = (X1+X2)/2
			Theta_Sgnl = XNew
			CALL FIND_OPT_PMF_1D(Theta_Idlr_Guess,FNew,Info)
			FNew = FNew - TargetPMF
			Error = DABS(X2-X1)
			IF (FNew .GT. 0d0) THEN
				Done = .TRUE.
			ELSE
				X2 = XNew
				F2 = FNew
			ENDIF
			Error = DABS(X2-X1)
			Iterations = Iterations + 1
		ENDDO

!		Step Through the interval looking for sign changes.  We'll keep
!		the answer closest to the original Theta_Sgnl
		XStep = (X1 - X2) / 12.00000001d0
		XStart = X2 + XStep
		XStop = X1
		FIn = F2
		DO XTest = XStart, XStop, XStep
			FOut = FIn
			Theta_Sgnl = XTest
			CALL FIND_OPT_PMF_1D(Theta_Idlr_Guess,FIn,Info)
			FIn = FIn - TargetPMF
			IF (DSIGN(1d0,FIn) .NE. DSIGN(1d0,FOut)) THEN
				X1 = XTest
				F1 = FIn
				X2 = XTest - XStep
				F2 = FOut
				IF (X1 .LT. X2) THEN
					XLow = X1
					XHigh = X2
				ELSE
					XLow = X2
					XHigh = X1
				ENDIF
				Error = 1d0
				Iterations = 0
!				Use Secant Method to find the value
				DO WHILE ((Error .GT. 1d-8) .AND. (Iterations .LT. 400))
					
					Theta_Sgnl = X2
					CALL FIND_OPT_PMF_1D(Theta_Idlr_Guess,F2,Info)
					F2 = F2 - TargetPMF

					XNew = X2 - F2*( X2-X1 ) / ( F2-F1 )

					IF (XNew .GT. XHigh) THEN
						XNew = XHigh
						Error = 1d0
					ELSEIF (XNew .LT. XLow) THEN
						XNew = XLow
						Error = 1d0
					ELSE
						Error = DABS(F2)
					ENDIF

					X1=X2
					F1 = F2
					X2 = XNew
					Iterations = Iterations + 1
				ENDDO
				IF (Iterations .NE. 400) THEN
					XAnswer = XNew
				ENDIF
			ENDIF
		ENDDO

		Theta_Sgnl = Theta_Sgnl_In
		Theta_Idlr = Theta_Idlr_In
	ENDIF
	Theta_AT_PMFVALUE = XAnswer

END FUNCTION THETA_AT_PMFVALUE

!***************************************************************************
!
!      FUNCTION PHI_AT_PMFVALUE(Direction, TargetPMF)
!
!		Looks for a Phi_Sgnl where the phase matching function 
!		matches a specified value.  Starts looking in the direction (+,-)
!		specified by the argument Direction.
!
!***************************************************************************
REAL*8 FUNCTION PHI_AT_PMFVALUE(Direction, TargetPMF)
	USE PM_DATA
	IMPLICIT NONE

!	Arguments

	REAL*8, INTENT(IN) :: Direction
	REAL*8, INTENT(IN) :: TargetPMF

!	Local Variables

	REAL*8 X1, X2, F1, F2, Temp, PHASE_MATCH_FUNC, XLow, XHigh
	REAL*8 Error, FNew, XNew, Phi_Sgnl_In, Max_Phi_Out
	INTEGER Info, Iterations
	LOGICAL Success, InRange

!	Begin Function
	Phi_Sgnl_In = Phi_Sgnl

	X1 = Phi_Sgnl
	X2 = Phi_Sgnl + DSIGN(Pi/2d0,Direction)
	F1 = PHASE_MATCH_FUNC(Success) - TargetPMF
	Phi_Sgnl = X2
	F2 = PHASE_MATCH_FUNC(Success) - TargetPMF

!	Make sure we have a point where the PMF Value is bigger than target
	IF (F1 .LT. 0d0) THEN
		PHI_AT_PMFVALUE = Phi_Sgnl
	ELSEIF (F2 .GT. 0d0) THEN
		PHI_AT_PMFVALUE = X2
	ELSE

!		Do a Binary Search to Bracket the point
		XNew = X2
		FNew = -1d0
		Iterations = 0
		DO WHILE ((FNew .LT. 0d0) .AND. (Iterations .LT. 400))
			X2 = XNew
			F2 = FNew
			XNew = (X1+X2)/2
			Phi_Sgnl = XNew
			FNew = PHASE_MATCH_FUNC(Success) - TargetPMF
			Iterations = Iterations + 1
		ENDDO
		X1 = XNew
		F1 = FNew


!		Use Secant Method to find the value
		IF (X1 .GT. X2) THEN
			XLow = X2
			XHigh = X1
		ELSE
			XLow = X1
			XHigh = X2
		ENDIF
		Error = 1d0
		Iterations = 0
		DO WHILE ((Error .GT. 1d-8) .AND. (Iterations .LT. 400))
			
			Phi_Sgnl = X2
			F2 = PHASE_MATCH_FUNC(Success) - TargetPMF

			XNew = X2 - F2*( X2-X1 ) / ( F2-F1 )

			IF (XNew .GT. XHigh) THEN
				XNew = XHigh
				Error = 1d0
			ELSEIF (XNew .LT. XLow) THEN
				XNew = XLow
				Error = 1d0
			ELSE
				Error = DABS(F2)
			ENDIF

			X1=X2
			F1 = F2
			X2 = XNew
			Iterations = Iterations + 1
		ENDDO

		PHI_AT_PMFVALUE = XNew
	ENDIF
	Phi_Sgnl = Phi_Sgnl_In
END FUNCTION PHI_AT_PMFVALUE

!***************************************************************************
!
!      FUNCTION THETA_IDLR_AT_PMFVALUE_(Direction, TargetPMF)
!
!		Looks for a Theta_Idlr where the phase matching function 
!		matches a specified value.  Starts looking in the direction (+,-)
!		specified by the argument Direction.
!
!***************************************************************************
REAL*8 FUNCTION THETA_IDLR_AT_PMFVALUE(Direction, TargetPMF)
	USE PM_DATA
	IMPLICIT NONE

!	Arguments

	REAL*8, INTENT(IN) :: Direction
	REAL*8, INTENT(IN) :: TargetPMF

!	Local Variables

	REAL*8 X1, X2, F1, F2, Temp, PHASE_MATCH_FUNC, XLow, XHigh
	REAL*8 Error, FNew, XNew, Theta_Idlr_In, Max_Theta_Out
	INTEGER Info, Iterations
	LOGICAL Success, InRange

!	Begin Function

	X1 = Theta_Idlr
	X2 = Theta_Idlr + DSIGN(Pi/2d0,Direction)

!	Make sure we have a point where the PMF Value is bigger than target
	F1 = PHASE_MATCH_FUNC(Success) - TargetPMF
	IF (F1 .GT. 0d0) THEN
		Theta_Idlr_In = Theta_Idlr

!		Do a Binary Search to Bracket the point
		XNew = X2
		FNew = -1d0
		Iterations = 0
		DO WHILE ((FNew .LT. 0d0) .AND. (Iterations .LT. 400))
			X2 = XNew
			F2 = FNew
			XNew = (X1+X2)/2
			Theta_Idlr = XNew
			FNew = PHASE_MATCH_FUNC(Success) - TargetPMF
			Iterations = Iterations + 1
		ENDDO
		X1 = XNew
		F1 = FNew


!		Use Secant Method to find the value
		IF (X1 .GT. X2) THEN
			XLow = X2
			XHigh = X1
		ELSE
			XLow = X1
			XHigh = X2
		ENDIF
		Error = 1d0
		Iterations = 0
		DO WHILE ((Error .GT. 1d-8) .AND. (Iterations .LT. 400))
			
			Theta_Idlr = X2
			F2 = PHASE_MATCH_FUNC(Success) - TargetPMF

			XNew = X2 - F2*( X2-X1 ) / ( F2-F1 )

			IF (XNew .GT. XHigh) THEN
				XNew = XHigh
				Error = 1d0
			ELSEIF (XNew .LT. XLow) THEN
				XNew = XLow
				Error = 1d0
			ELSE
				Error = DABS(F2)
			ENDIF

			X1=X2
			F1 = F2
			X2 = XNew
			Iterations = Iterations + 1
		ENDDO

		Theta_Idlr = Theta_Idlr_In
		Theta_IDLR_AT_PMFVALUE = XNew
	ELSE
		Theta_IDLR_AT_PMFVALUE = Theta_Idlr
	ENDIF
END FUNCTION THETA_IDLR_AT_PMFVALUE


!***************************************************************************
!
!      FUNCTION PMF_INTEGRAL(Success)
!
!***************************************************************************

REAL*8 FUNCTION PMF_INTEGRAL(Success)
	USE CRYSTAL_DATA; USE PM_DATA
	IMPLICIT NONE

!	Arguments

	LOGICAL, INTENT(INOUT) :: Success

!	Local Variables
	INTEGER Info, XIndex, YIndex, Resolution
	REAL*8 Delta_Phi, Delta_Theta, PMF_Value, Phi_Idlr_In, Theta_Idlr_In
	REAL*8 Phi_Idlr_Min, Phi_Idlr_Max, Phi_Idlr_Step, Integral
	REAL*8 Theta_Idlr_Min, Theta_Idlr_Max, Theta_Idlr_Step
	REAL*8 PHASE_MATCH_FUNC, PHI_AT_PMFVALUE, THETA_IDLR_AT_PMFVALUE

!	Begin Subroutine		
	Phi_Idlr_In = Phi_Idlr
	Theta_Idlr_In = Theta_Idlr

!	Decide on the limits
	CALL FIND_OPT_PMF_1D(Theta_Idlr_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
!	Phi_Idlr_Min = Phi_Idlr-Delta_Phi
!	Phi_Idlr_Max = 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
	Theta_Idlr_Min = Theta_Idlr-Delta_Theta
	Theta_Idlr_Max = Theta_Idlr+Delta_Theta

!	Perform the Integral
	Integral = 0d0
	Resolution = 100
!	Phi_Idlr_Step = (Phi_Idlr_Max - Phi_Idlr_Min) / (Resolution - 1)
	Theta_Idlr_Step = (Theta_Idlr_Max - Theta_Idlr_Min) / (Resolution - 1)

!	DO XIndex = 1, Resolution, 1
!		Phi_Idlr = Phi_Idlr_Min + (XIndex - 1)*Phi_Idlr_Step 
	
		DO YIndex = 1, Resolution, 1
			Theta_Idlr = Theta_Idlr_Min + (YIndex - 1)*Theta_Idlr_Step

			PMF_Value = PHASE_MATCH_FUNC(Success)

			Integral = Integral + PMF_Value*Theta_Idlr_Step !*DSIN(Theta_Idlr)*Phi_Idlr_Step
	
		END DO
!	END DO

	Phi_Idlr = Phi_Idlr_In
	Theta_Idlr = Theta_Idlr_In

	PMF_INTEGRAL = Integral
END FUNCTION PMF_INTEGRAL