!************************************************************************************** ! ! 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