PROGRAM PHASEMATCHMAIN IMPLICIT NONE INTEGER CRYSTAL,I INTEGER ITER INTEGER OUTCHAN,INCHAN INTEGER MENU,MENU3 INTEGER TYPEPM INTEGER N,LW,INFO PARAMETER (N=2) REAL*8 PI,W(N*(N+10)),F,X0(N),FACTOR,XINIT1,XINIT2 REAL*8 LAMBDA,LAMBDA2,LAMBDA_IDLER,LAMBDA_PUMP,LAMBDA_SIGNAL,LAMBDA_START,LAMBDA_STOP REAL*8 LAMBDA_INF,LAMBDA_MIN,LAMBDA_MAX,LAMBDA_MID1,LAMBDA_MID2 REAL*8 LAMBDA_DEG,LAMBDA_CONJ,LAMBDA_INC REAL*8 ABS_MIN,ABS_MAX REAL*8 N_PUMP,N_SIG,N_IDLER REAL*8 NF,NS,PS,QQ ! theta of signal (1) and idler (2) in pump coordinate system (x",y",z") REAL*8 ALPHA,BETA ! alpha is phi-pump and beta is theta-pump REAL*8 PHI_SIG2,PHI_SIGNAL,PHI_IDLER,PHI_INC REAL*8 THETA_INC REAL*8 THETA(2) REAL*8 THETA_START,THETA_STOP,PHI_START,PHI_STOP REAL*8 THETA_SIGNAL_START1,THETA_SIGNAL_STOP1,THETA_IDLER_START1,THETA_IDLER_STOP1 REAL*8 THETA_SIGNAL_START2,THETA_SIGNAL_STOP2,THETA_IDLER_START2,THETA_IDLER_STOP2 REAL*8 THETA_MIN,THETA_MAX,PHI_MIN,PHI_MAX REAL*8 DT_SIG_P, DT_IDL_P,DT_SIG_N,DT_IDL_N,DT_SIG_P_ext,DT_SIG_N_ext REAL*8 ANG1,ANG2,ANG1ext,ANG2ext,ANGSIG,ANG1MAX,ANG1MIN,ANGIDL,ANG1P,ANG1N,ANG1_INF REAL*8 TX,TY,TZ,NX,NY,NZ,TX_PUMP,TY_PUMP,TZ_PUMP REAL*8 XI(2,2) REAL*8 CryLength REAL*8 Pump_W,Pump_FWHM REAL*8 DELTA_K REAL*8 DKxyz(3) REAL*8 TEMP,TEMP2,LOOP1P,LOOP1N REAL*8 INC REAL*8 DOWN_INT,DOWN_INT_INC REAL*8 DK REAL*8 P(N) REAL*8 ITER1,ITER2 REAL*8 COUNTER REAL*8 DUMMY2 REAL*8 FTOL REAL*8 TOLERANCE PARAMETER (PI = 3.141592653589793238) PARAMETER (FTOL=1.E-10) PARAMETER (TOLERANCE=1.E-10) PARAMETER (XINIT1=0.03) PARAMETER (XINIT2=0.03) CHARACTER*30 DUMMY CHARACTER*30 PROGRAMNAME CHARACTER*30 FILENAME COMMON /CM/ CRYSTAL, TYPEPM, TX_PUMP,TY_PUMP,TZ_PUMP, & N_PUMP, N_SIG, N_IDLER, ALPHA, BETA, PHI_IDLER, PHI_SIGNAL, & LAMBDA_IDLER,LAMBDA_SIGNAL,LAMBDA_PUMP, DKxyz COMMON /CM2/ PROGRAMNAME EXTERNAL DeltaK ! CALL F_ABOUTSTRING ('Non-collinear Phase-matching','by Alan Migdall', 'with N. Boeuf, D. Branning, I. Chaperot, E. Dauler,','S. Guerin, G. Jaeger, A. Muller','NIST Optical Technology Division') PROGRAMNAME=" PHASEMATCH V 4 (08/20/99) " MENU = 0 LAMBDA=1 DO WHILE ( MENU .NE. 12 ) WRITE(*, *) PROGRAMNAME ! WRITE(*, *) "" WRITE(*, *) '(1) 2D Plot, Nx,Ny,Nz = f(Lambda)' ! WRITE(*, *) "" WRITE(*, *) '(2) 3D Plot, nSlow - nFast = f(Theta,Phi)' ! WRITE(*, *) "" WRITE(*, *) '(3) 3D Plot, minimum DeltaK = f(Theta,Phi)' WRITE(*, *) "" WRITE(*, *) '(4) 3D Plot, Phase-matching function = f(Dk-transverse,Dkz)' WRITE(*, *) "" WRITE(*, *) '(5) 2D Plot, Theta signal vs. Theta idler (lambda signal fixed)' WRITE(*, *) ' for a chosen value of the Phase-matching function' WRITE(*, *) "" WRITE(*, *) '(6) Polar plot, (Optimum Theta Idler,Optimum Theta Signal) = f(Phi)' WRITE(*, *) "" WRITE(*, *) '(7) 2D Plot, Optimum Theta = f(lambda) at chosen Phi' WRITE(*, *) ' with spreading in Theta,Phi Theta idler fixed' WRITE(*, *) "" WRITE(*, *) '(8) 2D Plot, Optimum Theta = f(lambda) at chosen Phi' WRITE(*, *) ' with spreading in Theta,Phi' WRITE(*, *) "" WRITE(*, *) '(9) 3D plot, Phase-matching function = f(Theta,Lambda)' WRITE(*, *) "" WRITE(*, *) '(10) 3D plot, Phase-matching function = f(Phi,Lambda)' ! WRITE(*, *) "" WRITE(*, *) '(11) HELP' ! WRITE(*, *) "" WRITE(*, *) '(12) Stop' ! WRITE(*, *) "" WRITE(*, "(' Choice:',$)") READ(*, *) MENU IF ((MENU .LE. 10).AND.(MENU .GE. 1)) THEN CALL INITFILE(FILENAME,OUTCHAN,INCHAN,CRYSTAL,MENU) IF ((MENU .NE. 4).AND.(MENU .NE. 1)) THEN WRITE(*, *) 'Enter scale factor to change the initial guess' READ(INCHAN,*) FACTOR ENDIF ENDIF IF ((MENU .LE. 10).AND.(MENU .GE. 7)) THEN MENU3 = 0 WRITE(*, *) 'How will the initial guess for the calculations be chosen?' WRITE(*, *) "" WRITE(*, *) '(1) Use same initial guess for each downconversion opening angle THETA.' WRITE(*, *) "" WRITE(*, *) '(2) Use previous result for theta optimum as new initial guess.' WRITE(*, *) "" WRITE(*, *) 'Choice:' READ(*, *) MENU3 ENDIF !*************************************************************************** IF ( MENU .EQ. 1 ) THEN ! CALL OutwindowClear WRITE(*, *) '(1) 2D Plot, Nx,Ny,Nz = f(Lambda)' WRITE(OUTCHAN, *) '(1) 2D Plot, Nx,Ny,Nz = f(Lambda)' WRITE(OUTCHAN, *) ' Lambda, Nx, Ny, Nz' CALL INDEX(CRYSTAL, LAMBDA, NX, NY, NZ, ABS_MIN, ABS_MAX) DO LAMBDA = ABS_MIN ,ABS_MAX , (ABS_MAX - ABS_MIN)/100. CALL INDEX(CRYSTAL,LAMBDA,NX,NY,NZ, ABS_MIN,ABS_MAX) WRITE(OUTCHAN,'(4(F10.5,'',''))') LAMBDA,NX,NY,NZ ENDDO CLOSE(INCHAN) CLOSE(OUTCHAN) pause !*************************************************************************** ELSE IF (MENU .EQ. 2) THEN ! CALL OutwindowClear WRITE(*, *) '(2) 3D Plot, nSlow - nFast = f(Theta,Phi)' WRITE(OUTCHAN, *) '(2) 3D Plot, nSlow - nFast = f(Theta,Phi)' WRITE(*, "(' Pump Wavelength (microns) =', $)") READ(INCHAN, *) LAMBDA WRITE(OUTCHAN,"(' Pump Wavelength (microns) =', F10.5)") LAMBDA WRITE(*,"(' Phi pump increment (deg) =',$)") READ(INCHAN,*) PHI_INC WRITE(OUTCHAN,"(' Phi pump increment (deg) =', F10.1)") PHI_INC PHI_INC = PHI_INC * PI/ 180. WRITE(*,"(' Theta pump increment (deg) =',$)") READ(INCHAN,*) THETA_INC WRITE(OUTCHAN,"(' Theta pump increment (deg) =', F10.1)") THETA_INC THETA_INC = THETA_INC * PI/ 180. WRITE(OUTCHAN,*) 'Phi pump , Theta pump , nSlow-nFast' ! We increment over the PHI (in the rotated coordinates system). DO ALPHA = 5*PI/180 , 2*PI- 5*PI/180 , PHI_INC WRITE( * , * ) "" WRITE( * , "(' Phi-pump =', F10.1)") ALPHA/PI*180. DO BETA = 5*PI/180 , PI- 5*PI/180 , THETA_INC ! The result is NaN when the angles are too close to 0 deg and 180 deg TX = COS(ALPHA)*SIN(BETA) TY = SIN(ALPHA)*SIN(BETA) TZ = COS(BETA) CALL INDEX(CRYSTAL, LAMBDA, NX,NY, NZ, ABS_MIN, ABS_MAX) PS = (TX**2)*(1/NY**2 + 1/NZ**2) + & (TY**2)*(1/NX**2 + 1/NZ**2) + & (TZ**2)*(1/NX**2 + 1/NY**2) QQ = (TX**2)/(NY**2 * NZ**2) + & (TY**2)/(NX**2 * NZ**2) + & (TZ**2)/(NX**2 * NY**2) NF = SQRT( 2.0/(PS + SQRT(PS**2 - 4*QQ)) ) NS = SQRT( 2.0/(PS - SQRT(PS**2 - 4*QQ)) ) WRITE(OUTCHAN,"(F10.1,',',F10.1,',',F10.7)") & ALPHA * 180./PI,BETA * 180./PI,NS - NF ENDDO ENDDO CLOSE(INCHAN) CLOSE(OUTCHAN) pause !*************************************************************************** ELSE IF ( MENU .EQ. 3 ) THEN ! CALL OutwindowClear WRITE(*, *) '(3) 3D Plot, minimum DeltaK = f(Theta,Phi)' WRITE(OUTCHAN, *) '(3) 3D Plot, minimum DeltaK = f(Theta,Phi)' WRITE(*, "(' Wavelength of Pump (microns) =', $)") READ(INCHAN, *) LAMBDA_PUMP WRITE(OUTCHAN,"(' Wavelength of Pump (microns) =', F10.5)") LAMBDA_PUMP WRITE(*, "(' Wavelength of Signal (microns) =', $)") READ(INCHAN, *) LAMBDA_SIGNAL WRITE(OUTCHAN,"(' Wavelength of Signal (microns) =', F10.5)") LAMBDA_SIGNAL LAMBDA_IDLER = 1/(1/LAMBDA_PUMP - 1/LAMBDA_SIGNAL) WRITE(OUTCHAN,"(' Wavelength of Idler (microns) =', F10.5)") LAMBDA_IDLER WRITE(*,"(' Wavelength of Idler (microns) =', F10.5)") LAMBDA_IDLER WRITE(OUTCHAN,*) ' Nx, Ny, Nz' CALL INDEX(CRYSTAL, LAMBDA_PUMP, NX,NY, NZ, ABS_MIN, ABS_MAX) WRITE(OUTCHAN,"('Pump: ',F10.5,F10.5,F10.5)") NX,NY,NZ CALL INDEX(CRYSTAL, LAMBDA_SIGNAL, NX,NY, NZ, ABS_MIN, ABS_MAX) WRITE(OUTCHAN,"('Signal:', F10.5,F10.5,F10.5)") NX,NY,NZ CALL INDEX(CRYSTAL, LAMBDA_IDLER, NX,NY, NZ, ABS_MIN, ABS_MAX) WRITE(OUTCHAN,"('Idler: ', F10.5,F10.5,F10.5)") NX,NY,NZ WRITE(*,"(' Phi pump increment (deg)=',$)") READ(INCHAN,*) PHI_INC WRITE(OUTCHAN,"(' Phi pump increment (deg)=', F10.1)") PHI_INC PHI_INC = PHI_INC * PI/ 180. WRITE(*,"(' Theta pump increment (deg)=',$)") READ(INCHAN,*) THETA_INC WRITE(OUTCHAN,"(' Theta pump increment (deg)=', F10.1)") THETA_INC THETA_INC = THETA_INC * PI/ 180. WRITE(*,*) 'Phase Matching TYPEPMs' WRITE(*,*) ' TYPEPM Pump Signal Idler' WRITE(*,*) ' 1 FAST SLOW SLOW' WRITE(*,*) ' 2 FAST FAST SLOW' WRITE(*,*) ' 3 FAST SLOW FAST' WRITE(*,"(' Choice: ',$)") READ(INCHAN,*) TYPEPM IF ( TYPEPM .EQ. 1 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Slow Idler=Slow' ELSE IF ( TYPEPM .EQ. 2 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Fast Idler=Slow' ELSE IF ( TYPEPM .EQ. 3 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Slow Idler=Fast' ENDIF WRITE(OUTCHAN, *) "PHI,THETA,delta-K" IF ( LAMBDA_SIGNAL .GE. ABS_MIN .AND. & LAMBDA_SIGNAL .LE. ABS_MAX .AND. & LAMBDA_PUMP .GE. ABS_MIN .AND. & LAMBDA_PUMP .LE. ABS_MAX .AND. & LAMBDA_IDLER .GE. ABS_MIN .AND. & LAMBDA_IDLER .LE. ABS_MAX ) THEN DO ALPHA = 5*PI/180 , 2*PI-5*PI/180, PHI_INC WRITE( * , * ) "" WRITE( * , "(' Phi-pump =',F10.1)") ALPHA/PI*180. DO BETA = 5*PI/180 , PI- 5*PI/180 , THETA_INC ! The result is NaN when the angles are too close to 0 deg and 180 deg CALL KPUMP LW=N*(N+10) ! Initial guesses for theta (signal,idler) optimum X0(1)=XINIT1*FACTOR X0(2)=XINIT2*FACTOR CALL UNCMND (N,X0,DeltaK,THETA,F,INFO,W,LW) WRITE(OUTCHAN,"(F10.1,',',F10.1,',',F10.7)") ALPHA * 180./PI,BETA * 180./PI,F WRITE( * , "(' Theta-pump =',F10.1 )") BETA/PI*180. ENDDO ENDDO ELSE WRITE( OUTCHAN , * ) " Not valid wavelength range " WRITE(*,*) " Not valid wavelength range " ENDIF ! if valid wavelength range WRITE(*, *) "Done with phase matching." CLOSE(INCHAN) CLOSE(OUTCHAN) pause !*************************************************************************** ELSE IF ( MENU .EQ. 4 ) THEN ! CALL OutwindowClear WRITE(*, *) '(4) 3D Plot, Phase-matching function = f(Dk-transverse,Dkz)' WRITE(OUTCHAN, *) '(4) 3D Plot, Phase-matching function = f(Dk-transverse,Dkz)' WRITE(*,"(' Length of Crystal (microns) =',$)") READ(INCHAN,*) CryLength WRITE(OUTCHAN,"(' Crystal length =',F10.0)") CryLength WRITE(*, "(' Gaussian Pump (intensity FWHM in microns) =', $)") READ(INCHAN, *) Pump_FWHM WRITE(OUTCHAN, "(' Gaussian Pump FWHM =', F10.0)") Pump_FWHM Pump_W=Pump_FWHM*0.8493218 WRITE(OUTCHAN,*) 'Dkz',',','Dk-transverse',',','Downconversion Intensity' DO TEMP=-3/Pump_W,3/Pump_W,3/Pump_W/50 WRITE(*,"(' Dkz =',F10.7 )") TEMP DO TEMP2=-3/Pump_W,3/Pump_W,3/Pump_W/50 DUMMY2=0.5*CryLength*TEMP DOWN_INT=EXP(-0.5*Pump_W**2*TEMP2**2)*(SIN(DUMMY2)/DUMMY2)**2 WRITE(OUTCHAN,'(3(F10.8,'',''))') TEMP,TEMP2,DOWN_INT ENDDO ENDDO CLOSE(INCHAN) CLOSE(OUTCHAN) pause !*************************************************************************** ELSE IF ( MENU .EQ. 5) THEN ! CALL OutwindowClear WRITE(*, *) '(5) 2D Plot, Theta signal vs. Theta idler (lambda signal fixed)' WRITE(*, *) ' for a chosen value of the Phase-matching function' WRITE(OUTCHAN, *) '(5) 2D Plot, Theta signal vs. Theta idler (lambda signal fixed)' WRITE(OUTCHAN, *) ' for a chosen value of the Phase-matching function' WRITE(*, "(' Wavelength of Pump (microns) =', $)") READ(INCHAN, *) LAMBDA_PUMP WRITE(OUTCHAN,"(' Wavelength of Pump (microns) =', F10.5)") LAMBDA_PUMP WRITE(*, "(' Wavelength of Signal (microns) =', $)") READ(INCHAN, *) LAMBDA2 WRITE(OUTCHAN, "(' Wavelength of Signal (microns) =', F10.5)") LAMBDA2 LAMBDA_IDLER = 1/(1/LAMBDA_PUMP - 1/LAMBDA2) WRITE(OUTCHAN,"(' Wavelength of Idler (microns) =', F10.5)") LAMBDA_IDLER WRITE(*,"(' Wavelength of Idler (microns) =', F10.5)") LAMBDA_IDLER WRITE(OUTCHAN,*) ' Nx, Ny, Nz' CALL INDEX(CRYSTAL, LAMBDA_PUMP, NX,NY, NZ, ABS_MIN, ABS_MAX) WRITE(OUTCHAN,"('Pump: ',F10.5,F10.5,F10.5)") NX,NY,NZ WRITE(*,"(' Phi Pump (degrees) =',$)") READ(INCHAN,*) ALPHA WRITE(*,"(' Theta Pump (degrees) =',$)") READ(INCHAN,*) BETA WRITE(OUTCHAN,"(' Phi Pump (deg)=',F10.2)") ALPHA WRITE(OUTCHAN,"(' Theta Pump (deg)=',F10.2)") BETA ALPHA = ALPHA * PI/ 180. BETA = BETA * PI/ 180. WRITE(*,"(' Phi Signal (degrees) =',$)") READ(INCHAN,*) PHI_SIG2 WRITE(OUTCHAN,"(' Phi Signal (deg)=',F10.2)") PHI_SIG2 PHI_SIG2 = PHI_SIG2 * PI /180 WRITE(*,*) 'Phase Matching TYPEPMs' WRITE(*,*) ' TYPEPM Pump Signal Idler' WRITE(*,*) ' 1 FAST SLOW SLOW' WRITE(*,*) ' 2 FAST FAST SLOW' WRITE(*,*) ' 3 FAST SLOW FAST' WRITE(*,"(' Choice: ',$)") READ(INCHAN,*) TYPEPM IF ( TYPEPM .EQ. 1 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Slow Idler=Slow' ELSE IF ( TYPEPM .EQ. 2 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Fast Idler=Slow' ELSE IF ( TYPEPM .EQ. 3 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Slow Idler=Fast' ENDIF WRITE(*,"(' Length of Crystal (microns) =',$)") READ(INCHAN,*) CryLength WRITE(OUTCHAN,"(' Crystal length =',F10.0)") CryLength WRITE(*, "(' Gaussian Pump (intensity FWHM in microns) =', $)") READ(INCHAN, *) Pump_FWHM WRITE(OUTCHAN, "(' Gaussian Pump FWHM =', F10.0)") Pump_FWHM Pump_W=Pump_FWHM*0.8493218 WRITE(*, "(' Downconversion Intensity =', $)") READ(INCHAN, *) DOWN_INT WRITE(OUTCHAN, "(' Downconversion Intensity =', F10.3)") DOWN_INT WRITE(*, "(' Theta increment (rd) =', $)") READ(INCHAN, *) INC WRITE(OUTCHAN, "(' Theta increment (rd) =', F10.8)") INC ANG1=XINIT1*FACTOR ANG2=XINIT2*FACTOR CALL KPUMP LAMBDA_SIGNAL=LAMBDA2 IF ( LAMBDA_SIGNAL .GE. ABS_MIN .AND. & LAMBDA_SIGNAL .LE. ABS_MAX .AND. & LAMBDA_PUMP .GE. ABS_MIN .AND. & LAMBDA_PUMP .LE. ABS_MAX .AND. & LAMBDA_IDLER .GE. ABS_MIN .AND. & LAMBDA_IDLER .LE. ABS_MAX ) THEN !***************************************************** ! CALCULATION ANGLES OF MAX INTENSITY PHI_SIGNAL = PHI_SIG2 PHI_IDLER = PHI_SIGNAL + PI LW=N*(N+10) ! Initial guesses for theta (signal,idler) optimum X0(1)=ANG1 X0(2)=ANG2 CALL UNCMND (N,X0,DeltaK,THETA,F,INFO,W,LW) ! Signal angle THETA(1)=ABS(THETA(1)) ANG1=THETA( 1 ) WRITE( * , "('Internal angle of max. signal intensity (deg, mrad) =',F10.4,F10.4)") ANG1*180 / PI,ANG1*1000 WRITE( OUTCHAN , "('Internal angle of max. signal intensity (deg, mrad) =',F10.4,F10.4)") ANG1*180 / PI,ANG1*1000 ! Idler angle THETA(2)=ABS(THETA(2)) ANG2=THETA( 2 ) WRITE( * , "('Internal angle of max. idler intensity (deg, mrad) =',F10.4,F10.4)") ANG2*180 / PI,ANG2*1000 WRITE( OUTCHAN , "('Internal angle of max. idler intensity (deg, mrad) =',F10.4,F10.4)") ANG2*180 / PI,ANG2*1000 !***************************************************** ! Angular Theta Spread with both sig. and idl. angles changing over entire range WRITE(OUTCHAN,*) "Delta Kz",",","Delta K-transverse",",","Theta signal",",","Theta idler" THETA(1)=ANG1 THETA(2)=ANG2 PHI_SIGNAL=PHI_SIG2 PHI_IDLER=PHI_SIGNAL+PI ! This loop will be finding ranges of phase-matching angles for all values of wavelength ! between the user's input value of lambda_signal and its conjugate, lambda_idler, ! provided that neither of these wavelengths fall outside the range of validity ! of the Sellmeir coefficients for the chosen crystal. ! Angular Spread with both sig. and idl. angles changing for fixed lambda THETA_START=PI THETA_STOP=PI THETA_MIN=PI THETA_MAX=0. ITER1=ANG2 ANG1P=ANG1 ANG1N=ANG1 DO WHILE ((abs((THETA_STOP-ANG1P)).GT.TOLERANCE).OR.((ABS(ANG1P-THETA_START)).GT.TOLERANCE)) THETA(2)=ITER1 THETA_SIGNAL_START1=ANG1P THETA_SIGNAL_STOP1=PI/2 DO WHILE ((THETA_SIGNAL_STOP1-THETA_SIGNAL_START1).GT.TOLERANCE) TEMP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2. THETA(1) = TEMP CALL DeltaK(N,THETA,DELTA_K) DUMMY2=0.5*CryLength*DKxyz(3) IF ((EXP(-0.5*Pump_W**2*(DKxyz(1)**2+DKxyz(2)**2))*(SIN(DUMMY2)/DUMMY2)**2).LT.DOWN_INT) THEN THETA_SIGNAL_STOP1=TEMP ELSE THETA_SIGNAL_START1=TEMP ENDIF ENDDO THETA_STOP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2 THETA_MAX=MAX(THETA_MAX,THETA_STOP) IF ((THETA_STOP-ANG1P).GT.TOLERANCE) THEN WRITE(OUTCHAN,'(5(F10.8,'',''))') DKxyz(3),(DKxyz(2)**2+DKxyz(1)**2)**.5,THETA_STOP*180/PI,THETA(2)*180/PI ENDIF THETA_SIGNAL_START1=0 THETA_SIGNAL_STOP1=ANG1P DO WHILE ((THETA_SIGNAL_STOP1-THETA_SIGNAL_START1).GT.TOLERANCE) TEMP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2. THETA(1) = TEMP CALL DeltaK(N,THETA,DELTA_K) DUMMY2=0.5*CryLength*DKxyz(3) IF ((EXP(-0.5*Pump_W**2*(DKxyz(1)**2+DKxyz(2)**2))*(SIN(DUMMY2)/DUMMY2)**2).LT.DOWN_INT) THEN THETA_SIGNAL_START1=TEMP ELSE THETA_SIGNAL_STOP1=TEMP ENDIF ENDDO THETA_START=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2 THETA_MIN=MIN(THETA_MIN,THETA_START) IF ((ANG1P-THETA_START).GT.TOLERANCE) THEN WRITE(OUTCHAN,'(5(F10.8,'',''))') DKxyz(3),(DKxyz(2)**2+DKxyz(1)**2)**.5,THETA_START*180/PI,THETA(2)*180/PI ENDIF ITER1=ITER1+INC ANG1P=(THETA_START+THETA_STOP)/2 ENDDO THETA_START=PI THETA_STOP=PI ITER1=ANG2 DO WHILE ((abs((THETA_STOP-ANG1N)).GT.TOLERANCE).OR.((ABS(ANG1N-THETA_START)).GT.TOLERANCE)) THETA(2)=ITER1 THETA_SIGNAL_START1=ANG1N THETA_SIGNAL_STOP1=PI/2 DO WHILE ((THETA_SIGNAL_STOP1-THETA_SIGNAL_START1).GT.TOLERANCE) TEMP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2. THETA(1) = TEMP CALL DeltaK(N,THETA,DELTA_K) DUMMY2=0.5*CryLength*DKxyz(3) IF ((EXP(-0.5*Pump_W**2*(DKxyz(1)**2+DKxyz(2)**2))*(SIN(DUMMY2)/DUMMY2)**2).LT.DOWN_INT) THEN THETA_SIGNAL_STOP1=TEMP ELSE THETA_SIGNAL_START1=TEMP ENDIF ENDDO THETA_STOP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2 THETA_MAX=MAX(THETA_MAX,THETA_STOP) IF ((THETA_STOP-ANG1N).GT.TOLERANCE) THEN WRITE(OUTCHAN,'(5(F10.8,'',''))') DKxyz(3),(DKxyz(2)**2+DKxyz(1)**2)**.5,THETA_STOP*180/PI,THETA(2)*180/PI ENDIF THETA_SIGNAL_START1=0 THETA_SIGNAL_STOP1=ANG1N DO WHILE ((THETA_SIGNAL_STOP1-THETA_SIGNAL_START1).GT.TOLERANCE) TEMP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2. THETA(1) = TEMP CALL DeltaK(N,THETA,DELTA_K) DUMMY2=0.5*CryLength*DKxyz(3) IF ((EXP(-0.5*Pump_W**2*(DKxyz(1)**2+DKxyz(2)**2))*(SIN(DUMMY2)/DUMMY2)**2).LT.DOWN_INT) THEN THETA_SIGNAL_START1=TEMP ELSE THETA_SIGNAL_STOP1=TEMP ENDIF ENDDO THETA_START=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2 THETA_MIN=MIN(THETA_MIN,THETA_START) IF ((ANG1N-THETA_START).GT.TOLERANCE) THEN WRITE(OUTCHAN,'(5(F10.8,'',''))') DKxyz(3),(DKxyz(2)**2+DKxyz(1)**2)**.5,THETA_START*180/PI,THETA(2)*180/PI ENDIF ITER1=ITER1-INC ANG1N=(THETA_START+THETA_STOP)/2 ENDDO ELSE WRITE(*,*) 'One of the wavelengths is out of the validity range' WRITE(OUTCHAN,*) 'One of the wavelengths is out of the validity range' ENDIF CLOSE(INCHAN) CLOSE(OUTCHAN) pause !*************************************************************************** ELSE IF ( MENU .EQ. 6 ) THEN ! CALL OutwindowClear WRITE(*, *) '(6) Polar plot, (Optimum Theta Idler,Optimum Theta Signal) = f(Phi)' WRITE(OUTCHAN, *) '(6) Polar plot, (Optimum Theta Idler,Optimum Theta Signal) = f(Phi)' WRITE(*, "(' Wavelength of Pump (microns) =', $)") READ(INCHAN, *) LAMBDA_PUMP WRITE(OUTCHAN,"(' Wavelength of Pump (microns) =', F10.5)") LAMBDA_PUMP WRITE(*, "(' Wavelength of Signal (microns) =', $)") READ(INCHAN, *) LAMBDA_SIGNAL WRITE(OUTCHAN,"(' Wavelength of Signal (microns) =', F10.5)") LAMBDA_SIGNAL LAMBDA_IDLER = 1/(1/LAMBDA_PUMP - 1/LAMBDA_SIGNAL) WRITE(OUTCHAN,"(' Wavelength of Idler (microns) =', F10.5)") LAMBDA_IDLER WRITE(*,"(' Wavelength of Idler (microns) =', F10.5)") LAMBDA_IDLER WRITE(OUTCHAN,*) ' Nx, Ny, Nz' CALL INDEX(CRYSTAL, LAMBDA_PUMP, NX,NY, NZ, ABS_MIN, ABS_MAX) WRITE(OUTCHAN,"('Pump: ',F10.5,F10.5,F10.5)") NX,NY,NZ CALL INDEX(CRYSTAL, LAMBDA_SIGNAL, NX,NY, NZ, ABS_MIN, ABS_MAX) WRITE(OUTCHAN,"('Signal:', F10.5,F10.5,F10.5)") NX,NY,NZ CALL INDEX(CRYSTAL, LAMBDA_IDLER, NX,NY, NZ, ABS_MIN, ABS_MAX) WRITE(OUTCHAN,"('Idler: ', F10.5,F10.5,F10.5)") NX,NY,NZ WRITE(*,"(' Phi Pump (deg) =',$)") READ(INCHAN,*) ALPHA WRITE(OUTCHAN,"(' Phi Pump (deg) =',F10.2)") ALPHA ALPHA = ALPHA * PI/ 180. WRITE(*,"(' Theta Pump (deg)=',$)") READ(INCHAN,*) BETA WRITE(OUTCHAN,"(' Theta Pump (deg) =',F10.2)") BETA BETA = BETA * PI/ 180. WRITE(*,"(' Phi increment (deg)=',$)") READ(INCHAN,*) PHI_INC WRITE(OUTCHAN,"(' Phi increment (deg) =', F10.2)") PHI_INC PHI_INC = PHI_INC * PI/ 180. WRITE(*,*) 'Phase Matching TYPEPMs' WRITE(*,*) ' TYPEPM Pump Signal Idler' WRITE(*,*) ' 1 FAST SLOW SLOW' WRITE(*,*) ' 2 FAST FAST SLOW' WRITE(*,*) ' 3 FAST SLOW FAST' WRITE(*,"(' Choice: ',$)") READ(INCHAN,*) TYPEPM IF ( TYPEPM .EQ. 1 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Slow Idler=Slow' ELSE IF ( TYPEPM .EQ. 2 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Fast Idler=Slow' ELSE IF ( TYPEPM .EQ. 3 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Slow Idler=Fast' ENDIF CALL INDEX(CRYSTAL, LAMBDA, NX,NY, NZ, ABS_MIN, ABS_MAX) IF ( LAMBDA_SIGNAL .GE. ABS_MIN .AND. & LAMBDA_SIGNAL .LE. ABS_MAX .AND. & LAMBDA_PUMP .GE. ABS_MIN .AND. & LAMBDA_PUMP .LE. ABS_MAX .AND. & LAMBDA_IDLER .GE. ABS_MIN .AND. & LAMBDA_IDLER .LE. ABS_MAX ) THEN IF ( TYPEPM .EQ. 1 .OR. TYPEPM .EQ. 2 .OR. TYPEPM .EQ.3 ) THEN CALL KPUMP WRITE(OUTCHAN,"('Tx pump =',F10.5)") TX_PUMP WRITE(OUTCHAN,"('Ty pump =',F10.5)") TY_PUMP WRITE(OUTCHAN,"('Tz pump =',F10.5)") TZ_PUMP WRITE(OUTCHAN,"('Phi signal,Theta sig internal,Theta sig external,',$)") WRITE(OUTCHAN,"('Phi idler,Theta idl internal,Theta idl external')") WRITE(*,*) "Phi, Delta K" ! Thisloop allows two different initial guesses (one small & one big) for each PHI ! This is needed because each PHI may have two legal THETA results DO I=0,1 DO COUNTER = 0.0, 2*PI, PHI_INC PHI_SIGNAL=COUNTER PHI_IDLER = PHI_SIGNAL + PI LW=N*(N+10) X0(1)=XINIT1*(FACTOR+I*9)/3 X0(2)=XINIT2*(FACTOR+I*9)/3 CALL UNCMND (N,X0,DeltaK,THETA,F,INFO,W,LW) IF ( THETA( 1 ) .LE. 0 ) THEN PHI_SIGNAL = PHI_SIGNAL - PI THETA( 1 ) = ABS(THETA( 1 )) ENDIF IF ( THETA( 2 ) .LE. 0 ) THEN PHI_IDLER = PHI_IDLER - PI THETA( 2 ) = ABS(THETA( 2 )) ENDIF WRITE(*,"(F10.2,F10.5)") PHI_SIGNAL*180/PI,F IF ((F) .LE.0.000001 ) THEN WRITE(OUTCHAN,"(F10.2,',',F10.7,',',F10.7,','F10.2,',',F10.7,',',F10.7)") PHI_SIGNAL*180./PI, & THETA(1)*180./PI,asin(sin(THETA(1))*N_SIG)*180/PI, & PHI_IDLER*180./PI,THETA(2)*180./PI,asin(sin(THETA(2))*N_SIG)*180/PI ENDIF ENDDO ENDDO ELSE CALL ERR_MSG(" No such TYPEPM.", 17) ENDIF ! TYPEPM. ELSE WRITE(*,*) 'One of the wavelengths is out of the validity range' WRITE(OUTCHAN,*) 'One of the wavelengths is out of the validity range' pause ENDIF ! if valid wavelength range WRITE(*, *) "Done with phase matching." CLOSE(INCHAN) CLOSE(OUTCHAN) pause !****************************************************************************** ELSE IF ( MENU .EQ. 7) THEN ! CALL OutwindowClear WRITE(*, *) '(7) 2D Plot, Optimum Theta = f(lambda) at chosen Phi' WRITE(*, *) ' with spreading in Theta,Phi Theta idler fixed' WRITE(OUTCHAN, *) '(7) 2D Plot, Optimum Theta = f(lambda) at chosen Phi' WRITE(OUTCHAN, *) ' with spreading in Theta,Phi Theta idler fixed' WRITE(*, "(' Wavelength of Pump (microns) =', $)") READ(INCHAN, *) LAMBDA_PUMP WRITE(OUTCHAN,"(' Wavelength of Pump (microns) =', F10.5)") LAMBDA_PUMP WRITE(OUTCHAN,*) ' Nx, Ny, Nz' CALL INDEX(CRYSTAL, LAMBDA_PUMP, NX,NY, NZ, ABS_MIN, ABS_MAX) WRITE(OUTCHAN,"('Pump: ',F10.5,F10.5,F10.5)") NX,NY,NZ WRITE(*,"(' Phi Pump (degrees) =',$)") READ(INCHAN,*) ALPHA WRITE(*,"(' Theta Pump (degrees) =',$)") READ(INCHAN,*) BETA WRITE(OUTCHAN,"(' Phi Pump (deg)=',F10.5)") ALPHA WRITE(OUTCHAN,"(' Theta Pump (deg)=',F10.5)") BETA ALPHA = ALPHA * PI/ 180. BETA = BETA * PI/ 180. WRITE(*,"(' Phi Signal (degrees) =',$)") READ(INCHAN,*) PHI_SIG2 WRITE(OUTCHAN,"(' Phi Signal (deg)=',F10.5)") PHI_SIG2 PHI_SIG2 = PHI_SIG2 * PI /180 WRITE(*,*) 'Phase Matching TYPEPMs' WRITE(*,*) ' TYPEPM Pump Signal Idler' WRITE(*,*) ' 1 FAST SLOW SLOW' WRITE(*,*) ' 2 FAST FAST SLOW' WRITE(*,*) ' 3 FAST SLOW FAST' WRITE(*,"(' Choice: ',$)") READ(INCHAN,*) TYPEPM IF ( TYPEPM .EQ. 1 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Slow Idler=Slow' ELSE IF ( TYPEPM .EQ. 2 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Fast Idler=Slow' ELSE IF ( TYPEPM .EQ. 3 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Slow Idler=Fast' ENDIF WRITE(*,"(' Length of Crystal (microns) =',$)") ! Crystal Length is assigned to the variable CryLength READ(INCHAN,*) CryLength WRITE(OUTCHAN,"(' Crystal length =',F10.0)") CryLength ! Gaussian pump FWHM assigned to Pump_W WRITE(*, "(' Gaussian Pump (intensity FWHM in microns) =', $)") READ(INCHAN, *) Pump_FWHM WRITE(OUTCHAN, "(' Gaussian Pump FWHM =', F10.0)") Pump_FWHM ! See #6, p. 47 (Pump_W is now the amplitude sigma, rather than the intensity FWHM) Pump_W=Pump_FWHM*0.8493218 WRITE(*, "(' Lambda Start (microns) =', $)") READ(INCHAN, *) LAMBDA_START WRITE(OUTCHAN, "(' Lambda Start (microns) =', F10.5)") LAMBDA_START WRITE(*, "(' Lambda Stop (microns) =', $)") READ(INCHAN, *) LAMBDA_STOP WRITE(OUTCHAN, "(' Lambda Stop (microns) =', F10.5)") LAMBDA_STOP CALL KPUMP WRITE(*, "(' Downconversion Intensity =', $)") READ(INCHAN, *) DOWN_INT WRITE(OUTCHAN, "(' Downconversion Intensity =', F10.3)") DOWN_INT !***************************************************** ! Angular Spread with both sig. and idl. angles changing over entire range of Abs_min to Abs_max WRITE(OUTCHAN, *) "All angles following are in degrees." WRITE(OUTCHAN, "('LAMBDA SIGNAL,Theta Signal internal,internal positive spread,internal negative spread,',$)") WRITE(OUTCHAN, "('Theta Signal external,external positive spread,external negative spread,',$)") WRITE(OUTCHAN, "('Phi positive spread,Phi negative spread')") PHI_SIGNAL=PHI_SIG2 PHI_IDLER=PHI_SIGNAL+PI ! This loop will be finding ranges of phase-matching angles for all values of wavelength ! between the user's input value of lambda_signal and its conjugate, lambda_idler, ! provided that neither of these wavelengths fall outside the range of validity ! of the Sellmeir coefficients for the chosen crystal. WRITE(*,*) 'Lambda Signal','Theta signal opt','Delta K' ANG1=XINIT1*FACTOR ANG2=XINIT2*FACTOR F=1 DO LAMBDA_SIGNAL = LAMBDA_START, LAMBDA_STOP, (LAMBDA_STOP-LAMBDA_START)/500 LAMBDA_IDLER = 1/(1/LAMBDA_PUMP - 1/LAMBDA_SIGNAL) IF ( LAMBDA_SIGNAL .GE. ABS_MIN .AND. & LAMBDA_SIGNAL .LE. ABS_MAX .AND. & LAMBDA_PUMP .GE. ABS_MIN .AND. & LAMBDA_PUMP .LE. ABS_MAX .AND. & LAMBDA_IDLER .GE. ABS_MIN .AND. & LAMBDA_IDLER .LE. ABS_MAX ) THEN !***************************************************** ! CALCULATION OF THETA SIGNAL OPTIMUM AND THETA IDLER OPTIMUM ! Angular Spread with both sig. and idl. angles changing for fixed lambda ! set best guesses for optimum phase-matching angle at current lambdas; POWELL will find ! exact phase matching angles and return in THETA() LW=N*(N+10) IF ((MENU3==2).AND.(F<0.00001)) THEN X0(1)=ANG1 X0(2)=ANG2 ELSE X0(1)=XINIT1*FACTOR X0(2)=XINIT2*FACTOR ENDIF CALL UNCMND (N,X0,DeltaK,THETA,F,INFO,W,LW) ! Optimum Signal angle (rad) for starting point of iterations below ! THETA(1)=ABS(THETA(1)) ! makes sure all solutions found are positive theta values ANG1=THETA(1) ! ANG1 is Optimum Signal angle (rad): ! THETA(2)=ABS(THETA(2)) ANG2=THETA( 2 ) ! ANG2 is optimum idler angle (rad) IF (F<0.00001) THEN !***************************************************** ! CALCULATION OF THE SPREAD IN PHI SIGNAL WITH THETA SIGNAL=THETA OPTIMUM THETA(1)=ANG1 THETA(2)=ANG2 PHI_SIGNAL=PHI_SIG2 PHI_IDLER=PHI_SIG2+PI PHI_START=PHI_SIG2 PHI_STOP=PHI_SIG2+PI/2 DO WHILE ((PHI_STOP-PHI_START).GT.TOLERANCE) TEMP=(PHI_START+PHI_STOP)/2. PHI_SIGNAL = TEMP CALL DeltaK(N,THETA,DELTA_K) DUMMY2=0.5*CryLength*DKxyz(3) IF ((EXP(-0.5*Pump_W**2*(DKxyz(1)**2+DKxyz(2)**2))*(SIN(DUMMY2)/DUMMY2)**2).LT.DOWN_INT) THEN PHI_STOP=TEMP ELSE PHI_START=TEMP ENDIF ENDDO PHI_MAX=(PHI_START+PHI_STOP)/2 PHI_MIN=2*PHI_SIG2-PHI_MAX !***************************************************** ! CALCULATION OF THE SPREAD IN THETA SIGNAL THETA(1)=ANG1 THETA(2)=ANG2 PHI_SIGNAL=PHI_SIG2 PHI_IDLER=PHI_SIGNAL+PI THETA_SIGNAL_START2=0. THETA_SIGNAL_STOP2=1. ! intialize variables to execute loop the first time. ! This loop finds the highest value of the signal theta that can still produce a value of .5 for ! the phase-matching function, by allowing the idler theta value to vary. Each iteration uses the ! previous uppermost theta signal value to produce a new uppermost idler value, and then re-computes ! the uppermost theta signal value based on the new uppermost idler value. The loop terminates when ! the result converges to a value that is less than 1.E-12 radians away from the value obtained ! in the previous cycle. THETA_SIGNAL_START1=THETA(1) ! record the previous value of THETA(1), which resulted in a phase-matching function value of DOWN_INT (or 1 on 1st iteration) THETA_SIGNAL_STOP1=PI/2 ! the search will include all angles up to PI/2 THETA_SIGNAL_START2=THETA(1) ! this value will be compared to the new THETA_SIGNAL_STOP2 to test for DO WHILE (abs((THETA_SIGNAL_STOP1-THETA_SIGNAL_START1)).GT.TOLERANCE) TEMP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2. THETA(1) = TEMP CALL DeltaK(N,THETA,DELTA_K) DUMMY2=0.5*CryLength*DKxyz(3) IF ((EXP(-0.5*Pump_W**2*(DKxyz(1)**2+DKxyz(2)**2))*(SIN(DUMMY2)/DUMMY2)**2).LT.DOWN_INT) THEN THETA_SIGNAL_STOP1=TEMP ELSE THETA_SIGNAL_START1=TEMP ENDIF ENDDO IF (abs((THETA_SIGNAL_STOP1)).LT.TOLERANCE) THEN THETA(1)=0 ELSE THETA(1)=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2 ! after exiting loop, THETA(1) holds the value of the signal theta that results in phase-matching function = DOWN_INT ENDIF THETA_SIGNAL_STOP2=THETA(1) DT_SIG_P=THETA_SIGNAL_STOP2-ANG1 ! above lines compute the "error bars" or spreads, like HWHM, for plotting the ranges of allowed theta values ! This works in the same fashion os the positive spread, but the search is now conducted from 0 to the optimum ! theta values for signal and idler. THETA_SIGNAL_START2=1. THETA_SIGNAL_STOP2=0. THETA(1)=ANG1 THETA(2)=ANG2 THETA_SIGNAL_START1=-PI/2 THETA_SIGNAL_STOP1=THETA(1) THETA_SIGNAL_START2=THETA(1) DO WHILE (abs((THETA_SIGNAL_STOP1-THETA_SIGNAL_START1)).GT.TOLERANCE) TEMP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2. THETA(1) = TEMP CALL DeltaK(N,THETA,DELTA_K) DUMMY2=0.5*CryLength*DKxyz(3) IF ((EXP(-0.5*Pump_W**2*(DKxyz(1)**2+DKxyz(2)**2))*(SIN(DUMMY2)/DUMMY2)**2).LT.DOWN_INT) THEN THETA_SIGNAL_START1=TEMP ELSE THETA_SIGNAL_STOP1=TEMP ENDIF ENDDO IF (abs((THETA_SIGNAL_STOP1)).LT.TOLERANCE) THEN THETA(1)=0 ELSE THETA(1)=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2 ENDIF THETA_SIGNAL_STOP2=THETA(1) DT_SIG_N=ANG1-THETA_SIGNAL_STOP2 ! Conversion Internal angle -> External angle ! ANG1=abs(ANG1) ANG1ext = ( asin( sin(ANG1)*N_SIG ) ) DT_SIG_P_ext = ( ( asin( sin(ANG1ext+DT_SIG_P)*N_SIG ) ) - ( asin( sin(ANG1ext)*N_SIG ) ) ) DT_SIG_N_ext = ( ( asin( sin(ANG1ext)*N_SIG ) ) - ( asin( sin(ANG1ext-DT_SIG_N)*N_SIG ) ) ) WRITE(*,"(F10.5,F10.5,F10.5)") LAMBDA_SIGNAL,ANG1*180/PI,F WRITE(OUTCHAN,'(10(F20.7,'',''))') & LAMBDA_SIGNAL,ANG1*180/PI,DT_SIG_P*(180/PI),DT_SIG_N*(180/PI), & ANG1ext*180/PI,DT_SIG_P_ext*(180/PI),DT_SIG_N_ext*(180/PI), & (PHI_MAX-PHI_SIG2)*180/PI,(PHI_SIG2-PHI_MIN)*180/PI WRITE(*,*) '' ENDIF ENDIF ENDDO CLOSE(INCHAN) CLOSE(OUTCHAN) pause !*************************************************************************** ELSE IF ( MENU .EQ. 8 ) THEN ! CALL OutwindowClear WRITE(*, *) '(8) 2D Plot, Optimum Theta = f(lambda) at chosen Phi' WRITE(*, *) ' with spreading in Theta,Phi' WRITE(OUTCHAN, *) '(8) 2D Plot, Optimum Theta = f(lambda) at chosen Phi' WRITE(OUTCHAN, *) ' with spreading in Theta,Phi' WRITE(*, "(' Wavelength of Pump (microns) =', $)") READ(INCHAN, *) LAMBDA_PUMP WRITE(OUTCHAN,"(' Wavelength of Pump (microns) =', F10.5)") LAMBDA_PUMP WRITE(OUTCHAN,*) ' Nx, Ny, Nz' CALL INDEX(CRYSTAL, LAMBDA_PUMP, NX,NY, NZ, ABS_MIN, ABS_MAX) WRITE(OUTCHAN,"('Pump: ',F10.5,F10.5,F10.5)") NX,NY,NZ WRITE(*,"(' Phi Pump (degrees) =',$)") READ(INCHAN,*) ALPHA WRITE(*,"(' Theta Pump (degrees) =',$)") READ(INCHAN,*) BETA WRITE(OUTCHAN,"(' Phi Pump (deg)=',F10.2)") ALPHA WRITE(OUTCHAN,"(' Theta Pump (deg)=',F10.2)") BETA ALPHA = ALPHA * PI/ 180. BETA = BETA * PI/ 180. WRITE(*,"(' Phi Signal (degrees) =',$)") READ(INCHAN,*) PHI_SIG2 WRITE(OUTCHAN,"(' Phi Signal (deg)=',F10.2)") PHI_SIG2 PHI_SIG2 = PHI_SIG2 * PI /180 WRITE(*,*) 'Phase Matching TYPEPMs' WRITE(*,*) ' TYPEPM Pump Signal Idler' WRITE(*,*) ' 1 FAST SLOW SLOW' WRITE(*,*) ' 2 FAST FAST SLOW' WRITE(*,*) ' 3 FAST SLOW FAST' WRITE(*,"(' Choice: ',$)") READ(INCHAN,*) TYPEPM IF ( TYPEPM .EQ. 1 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Slow Idler=Slow' ELSE IF ( TYPEPM .EQ. 2 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Fast Idler=Slow' ELSE IF ( TYPEPM .EQ. 3 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Slow Idler=Fast' ENDIF WRITE(*,"(' Length of Crystal (microns) =',$)") READ(INCHAN,*) CryLength WRITE(OUTCHAN,"(' Crystal length =',F10.0)") CryLength WRITE(*, "(' Gaussian Pump (intensity FWHM in microns) =', $)") READ(INCHAN, *) Pump_FWHM WRITE(OUTCHAN, "(' Gaussian Pump FWHM =', F10.0)") Pump_FWHM Pump_W=Pump_FWHM*0.8493218 WRITE(*, "(' Lambda Start (microns) =', $)") READ(INCHAN, *) LAMBDA_START WRITE(OUTCHAN, "(' Lambda Start (microns) =', F10.5)") LAMBDA_START WRITE(*, "(' Lambda Stop (microns) =', $)") READ(INCHAN, *) LAMBDA_STOP WRITE(OUTCHAN, "(' Lambda Stop (microns) =', F10.5)") LAMBDA_STOP WRITE(*, "(' Downconversion Intensity =', $)") READ(INCHAN, *) DOWN_INT WRITE(OUTCHAN, "(' Downconversion Intensity =', F10.3)") DOWN_INT WRITE(*, "(' Theta increment (rd) =', $)") READ(INCHAN, *) INC WRITE(OUTCHAN, "(' Theta increment (rd) =', F10.8)") INC CALL KPUMP F=1 !***************************************************** ! Angular Spread with both sig. and idl. angles changing over entire range of Abs_min to Abs_max WRITE(OUTCHAN, *) "All angles following are in degrees." WRITE(OUTCHAN, "('LAMBDA SIGNAL,Theta Signal internal,internal positive spread,internal negative spread,',$)") WRITE(OUTCHAN, "('Theta Signal external,external positive spread,external negative spread,',$)") WRITE(OUTCHAN, "('Phi positive spread,Phi negative spread')") THETA(1)=ANG1 THETA(2)=ANG2 PHI_SIGNAL=PHI_SIG2 PHI_IDLER=PHI_SIGNAL+PI ! This loop will be finding ranges of phase-matching angles for all values of wavelength ! between the user's input value of lambda_signal and its conjugate, lambda_idler, ! provided that neither of these wavelengths fall outside the range of validity ! of the Sellmeir coefficients for the chosen crystal. ANG1=XINIT1*FACTOR ANG2=XINIT2*FACTOR WRITE(*,*) 'Lambda Signal','Theta signal opt','Delta K' DO LAMBDA_SIGNAL = LAMBDA_START, LAMBDA_STOP, (LAMBDA_STOP-LAMBDA_START)/200 LAMBDA_IDLER = 1/(1/LAMBDA_PUMP - 1/LAMBDA_SIGNAL) IF ( LAMBDA_SIGNAL .GE. ABS_MIN .AND. & LAMBDA_SIGNAL .LE. ABS_MAX .AND. & LAMBDA_PUMP .GE. ABS_MIN .AND. & LAMBDA_PUMP .LE. ABS_MAX .AND. & LAMBDA_IDLER .GE. ABS_MIN .AND. & LAMBDA_IDLER .LE. ABS_MAX ) THEN ! Angular Spread with both sig. and idl. angles changing for fixed lambda ! set best guesses for optimum phase-matching angle at current lambdas; POWELL will find ! exact phase matching angles and return in THETA() LW=N*(N+10) ! Initial guesses for theta (signal,idler) optimum IF ((MENU3==2).AND.(F<0.00001)) THEN X0(1)=ANG1 X0(2)=ANG2 ELSE X0(1)=XINIT1*FACTOR X0(2)=XINIT2*FACTOR ENDIF CALL UNCMND (N,X0,DeltaK,THETA,F,INFO,W,LW) ! Optimum Signal angle (rad) for starting point of iterations below ! THETA( 1 )=ABS(THETA( 1 )) ! makes sure all solutions found are positive theta values ANG1=THETA( 1 ) ! ANG1 is Optimum Signal angle (rad): ! THETA( 2 )=ABS(THETA( 2 )) ANG2=THETA( 2 ) ! ANG2 is optimum idler angle (rad) IF (F<0.00001) THEN !***************************************************** ! CALCULATION OF THE SPREAD IN PHI SIGNAL WITH THETA SIGNAL=THETA OPTIMUM THETA(1)=ANG1 THETA(2)=ANG2 PHI_SIGNAL=PHI_SIG2 PHI_IDLER=PHI_SIG2+PI PHI_START=PHI_SIG2 PHI_STOP=PHI_SIG2+PI/2 DO WHILE ((PHI_STOP-PHI_START).GT.TOLERANCE) THETA(1)=ANG1 THETA(2)=ANG2 TEMP=(PHI_START+PHI_STOP)/2. PHI_SIGNAL = TEMP CALL DeltaK(N,THETA,DELTA_K) DUMMY2=0.5*CryLength*DKxyz(3) IF ((EXP(-0.5*Pump_W**2*(DKxyz(1)**2+DKxyz(2)**2))*(SIN(DUMMY2)/DUMMY2)**2).LT.DOWN_INT) THEN PHI_STOP=TEMP ELSE PHI_START=TEMP ENDIF ENDDO PHI_MAX=(PHI_START+PHI_STOP)/2 PHI_MIN=2*PHI_SIG2-PHI_MAX !***************************************************** ! CALCULATION OF THE SPREAD IN THETA SIGNAL PHI_SIGNAL=PHI_SIG2 PHI_IDLER=PHI_SIG2+PI THETA(1)=ANG1 THETA(2)=ANG2 THETA_START=PI THETA_STOP=PI THETA_MIN=PI THETA_MAX=-PI ITER1=ANG2 ANG1P=ANG1 ANG1N=ANG1 DO WHILE ((abs((THETA_STOP-ANG1P)).GT.TOLERANCE).OR.((ABS(ANG1P-THETA_START)).GT.TOLERANCE)) THETA(2)=ITER1 THETA_SIGNAL_START1=ANG1P THETA_SIGNAL_STOP1=PI/2 DO WHILE (abs((THETA_SIGNAL_STOP1-THETA_SIGNAL_START1)).GT.TOLERANCE) TEMP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2. THETA(1) = TEMP CALL DeltaK(N,THETA,DELTA_K) DUMMY2=0.5*CryLength*DKxyz(3) IF ((EXP(-0.5*Pump_W**2*(DKxyz(1)**2+DKxyz(2)**2))*(SIN(DUMMY2)/DUMMY2)**2).LT.DOWN_INT) THEN THETA_SIGNAL_STOP1=TEMP ELSE THETA_SIGNAL_START1=TEMP ENDIF ENDDO THETA_STOP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2 THETA_MAX=MAX(THETA_MAX,THETA_STOP) THETA_SIGNAL_START1=-PI/2 THETA_SIGNAL_STOP1=ANG1P DO WHILE (abs((THETA_SIGNAL_STOP1-THETA_SIGNAL_START1)).GT.TOLERANCE) TEMP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2. THETA(1) = TEMP CALL DeltaK(N,THETA,DELTA_K) DUMMY2=0.5*CryLength*DKxyz(3) IF ((EXP(-0.5*Pump_W**2*(DKxyz(1)**2+DKxyz(2)**2))*(SIN(DUMMY2)/DUMMY2)**2).LT.DOWN_INT) THEN THETA_SIGNAL_START1=TEMP ELSE THETA_SIGNAL_STOP1=TEMP ENDIF ENDDO THETA_START=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2 THETA_MIN=MIN(THETA_MIN,THETA_START) ITER1=ITER1+INC ANG1P=(THETA_START+THETA_STOP)/2 ENDDO THETA_START=PI THETA_STOP=PI ITER1=ANG2 DO WHILE ((abs((THETA_STOP-ANG1N)).GT.TOLERANCE).OR.((ABS(ANG1N-THETA_START)).GT.TOLERANCE)) THETA(2)=ITER1 THETA_SIGNAL_START1=ANG1N THETA_SIGNAL_STOP1=PI/2 DO WHILE ((THETA_SIGNAL_STOP1-THETA_SIGNAL_START1).GT.TOLERANCE) TEMP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2. THETA(1) = TEMP CALL DeltaK(N,THETA,DELTA_K) DUMMY2=0.5*CryLength*DKxyz(3) IF ((EXP(-0.5*Pump_W**2*(DKxyz(1)**2+DKxyz(2)**2))*(SIN(DUMMY2)/DUMMY2)**2).LT.DOWN_INT) THEN THETA_SIGNAL_STOP1=TEMP ELSE THETA_SIGNAL_START1=TEMP ENDIF ENDDO THETA_STOP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2 THETA_MAX=MAX(THETA_MAX,THETA_STOP) THETA_SIGNAL_START1=-PI/2 THETA_SIGNAL_STOP1=ANG1N DO WHILE (abs((THETA_SIGNAL_STOP1-THETA_SIGNAL_START1)).GT.TOLERANCE) TEMP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2. THETA(1) = TEMP CALL DeltaK(N,THETA,DELTA_K) DUMMY2=0.5*CryLength*DKxyz(3) IF ((EXP(-0.5*Pump_W**2*(DKxyz(1)**2+DKxyz(2)**2))*(SIN(DUMMY2)/DUMMY2)**2).LT.DOWN_INT) THEN THETA_SIGNAL_START1=TEMP ELSE THETA_SIGNAL_STOP1=TEMP ENDIF ENDDO THETA_START=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2 THETA_MIN=MIN(THETA_MIN,THETA_START) ITER1=ITER1-INC ANG1N=(THETA_START+THETA_STOP)/2 ENDDO ! Conversion Internal angle -> External angle ! ANG1=abs(ANG1) ANG1ext = asin(sin(ANG1)*N_SIG) DT_SIG_P_ext = asin(sin(THETA_MAX)*N_SIG)-ANG1ext DT_SIG_N_ext = ANG1ext-asin(sin(THETA_MIN)*N_SIG) WRITE(*,"(F10.5,F10.5,F10.5)") LAMBDA_SIGNAL,ANG1*180/PI,F WRITE(OUTCHAN,'(10(F20.7,'',''))') & LAMBDA_SIGNAL,ANG1*180/PI,(THETA_MAX-ANG1)*180/PI,(ANG1-THETA_MIN)*180/PI, & ANG1ext*180/PI,DT_SIG_P_ext*180/PI,DT_SIG_N_ext*180/PI, & (PHI_MAX-PHI_SIG2)*180/PI,(PHI_SIG2-PHI_MIN)*180/PI ENDIF ENDIF ENDDO CLOSE(INCHAN) CLOSE(OUTCHAN) pause !*************************************************************************** ELSE IF ( MENU .EQ. 9 ) THEN ! CALL OutwindowClear WRITE(*,*) '(9) 3D plot, Phase-matching function = f(Theta,Lambda)' WRITE(OUTCHAN, *) '(9) 3D plot, Phase-matching function = f(Theta,Lambda)' WRITE(*, "(' Wavelength of Pump (microns) =', $)") READ(INCHAN, *) LAMBDA_PUMP WRITE(OUTCHAN,"(' Wavelength of Pump (microns) =', F10.5)") LAMBDA_PUMP WRITE(OUTCHAN,*) ' Nx, Ny, Nz' CALL INDEX(CRYSTAL, LAMBDA_PUMP, NX,NY, NZ, ABS_MIN, ABS_MAX) WRITE(OUTCHAN,"('Pump: ',F10.5,F10.5,F10.5)") NX,NY,NZ WRITE(*,"(' Phi Pump (degrees) =',$)") READ(INCHAN,*) ALPHA WRITE(OUTCHAN,"(' Phi Pump (deg)=',F10.5)") ALPHA WRITE(*,"(' Theta Pump (degrees) =',$)") READ(INCHAN,*) BETA WRITE(OUTCHAN,"(' Theta Pump (deg)=',F10.5)") BETA WRITE(*,"(' Phi Signal (degrees) =',$)") READ(INCHAN,*) PHI_SIG2 WRITE(OUTCHAN,"(' Phi Signal (deg)=',F10.5)") PHI_SIG2 ALPHA = ALPHA * PI/ 180. BETA = BETA * PI/ 180. PHI_SIG2= PHI_SIG2 * PI /180. WRITE(*,*) 'Phase Matching TYPEPMs' WRITE(*,*) ' TYPEPM Pump Signal Idler' WRITE(*,*) ' 1 FAST SLOW SLOW' WRITE(*,*) ' 2 FAST FAST SLOW' WRITE(*,*) ' 3 FAST SLOW FAST' WRITE(*,"(' Choice: ',$)") READ(INCHAN,*) TYPEPM IF ( TYPEPM .EQ. 1 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Slow Idler=Slow' ELSE IF ( TYPEPM .EQ. 2 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Fast Idler=Slow' ELSE IF ( TYPEPM .EQ. 3 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Slow Idler=Fast' ENDIF WRITE(*,"(' Length of Crystal (microns) =',$)") READ(INCHAN,*) CryLength WRITE(OUTCHAN,"(' Crystal length =',F10.0)") CryLength WRITE(*, "(' Gaussian Pump (intensity FWHM in microns) =', $)") READ(INCHAN, *) Pump_FWHM WRITE(OUTCHAN, "(' Gaussian Pump FWHM =', F10.0)") Pump_FWHM Pump_W=Pump_FWHM*0.8493218 WRITE(*, "(' Lambda Start (microns) =', $)") READ(INCHAN, *) LAMBDA_START WRITE(OUTCHAN, "(' Lambda Start (microns) =', F10.5)") LAMBDA_START WRITE(*, "(' Lambda Stop (microns) =', $)") READ(INCHAN, *) LAMBDA_STOP WRITE(OUTCHAN, "(' Lambda Stop (microns) =', F10.5)") LAMBDA_STOP WRITE(*, "(' Downconversion Intensity inc =', $)") READ(INCHAN, *) DOWN_INT_INC WRITE(OUTCHAN, "(' Downconversion Intensity inc =', F10.3)") DOWN_INT_INC WRITE(*, "(' Theta increment (rd) =', $)") READ(INCHAN, *) INC WRITE(OUTCHAN, "(' Theta increment (rd) =', F10.8)") INC WRITE(OUTCHAN, "('Lambda Signal,Theta Signal internal,Theta Signal external,',$)") WRITE(OUTCHAN, "('Downconversion Intensity')") ANG1=XINIT1*FACTOR ANG2=XINIT2*FACTOR CALL KPUMP WRITE(*,*) 'Lambda signal = ','Theta signal optimum = ' F=1 DO LAMBDA_SIGNAL = LAMBDA_START, LAMBDA_STOP, (LAMBDA_STOP-LAMBDA_START)/100 LAMBDA_IDLER=1/(1/LAMBDA_PUMP-1/LAMBDA_SIGNAL) IF ( LAMBDA_SIGNAL .GE. ABS_MIN .AND. & LAMBDA_SIGNAL .LE. ABS_MAX .AND. & LAMBDA_PUMP .GE. ABS_MIN .AND. & LAMBDA_PUMP .LE. ABS_MAX .AND. & LAMBDA_IDLER .GE. ABS_MIN .AND. & LAMBDA_IDLER .LE. ABS_MAX ) THEN PHI_SIGNAL=PHI_SIG2 PHI_IDLER=PHI_SIGNAL+PI LW=N*(N+10) ! Initial guesses for theta (signal,idler) optimum IF ((MENU3==2).AND.(F<0.00001)) THEN X0(1)=ANG1 X0(2)=ANG2 ELSE X0(1)=XINIT1*FACTOR X0(2)=XINIT2*FACTOR ENDIF CALL UNCMND (N,X0,DeltaK,THETA,F,INFO,W,LW) ! THETA( 1 )=ABS(THETA( 1 )) ANG1=THETA( 1 ) ! THETA( 2 )=ABS(THETA( 2 )) ANG2=THETA( 2 ) IF (F<0.00001) THEN WRITE(*,"(F10.5,F10.5)") LAMBDA_SIGNAL,THETA(1)*180/PI DO DOWN_INT=1,0.1,-DOWN_INT_INC THETA_START=PI THETA_STOP=PI THETA_MIN=PI THETA_MAX=-PI ITER1=ANG2 ANG1P=ANG1 ANG1N=ANG1 DO WHILE ((abs((THETA_STOP-ANG1P)).GT.TOLERANCE).OR.((ABS(ANG1P-THETA_START)).GT.TOLERANCE)) THETA(2)=ITER1 THETA_SIGNAL_START1=ANG1P THETA_SIGNAL_STOP1=PI/2 DO WHILE ((abs(THETA_SIGNAL_STOP1-THETA_SIGNAL_START1)).GT.TOLERANCE) TEMP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2. THETA(1) = TEMP CALL DeltaK(N,THETA,DELTA_K) DUMMY2=0.5*CryLength*DKxyz(3) IF ((EXP(-0.5*Pump_W**2*(DKxyz(1)**2+DKxyz(2)**2))*(SIN(DUMMY2)/DUMMY2)**2).LT.DOWN_INT) THEN THETA_SIGNAL_STOP1=TEMP ELSE THETA_SIGNAL_START1=TEMP ENDIF ENDDO THETA_STOP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2 THETA_MAX=MAX(THETA_MAX,THETA_STOP) THETA_SIGNAL_START1=-PI/2 THETA_SIGNAL_STOP1=ANG1P DO WHILE ((abs(THETA_SIGNAL_STOP1-THETA_SIGNAL_START1)).GT.TOLERANCE) TEMP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2. THETA(1) = TEMP CALL DeltaK(N,THETA,DELTA_K) DUMMY2=0.5*CryLength*DKxyz(3) IF ((EXP(-0.5*Pump_W**2*(DKxyz(1)**2+DKxyz(2)**2))*(SIN(DUMMY2)/DUMMY2)**2).LT.DOWN_INT) THEN THETA_SIGNAL_START1=TEMP ELSE THETA_SIGNAL_STOP1=TEMP ENDIF ENDDO THETA_START=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2 THETA_MIN=MIN(THETA_MIN,THETA_START) ITER1=ITER1+INC ANG1P=(THETA_START+THETA_STOP)/2 ENDDO THETA_START=PI THETA_STOP=PI ITER1=ANG2 DO WHILE ((abs((THETA_STOP-ANG1N)).GT.TOLERANCE).OR.((ABS(ANG1N-THETA_START)).GT.TOLERANCE)) THETA(2)=ITER1 THETA_SIGNAL_START1=ANG1N THETA_SIGNAL_STOP1=PI/2 DO WHILE ((abs(THETA_SIGNAL_STOP1-THETA_SIGNAL_START1)).GT.TOLERANCE) TEMP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2. THETA(1) = TEMP CALL DeltaK(N,THETA,DELTA_K) DUMMY2=0.5*CryLength*DKxyz(3) IF ((EXP(-0.5*Pump_W**2*(DKxyz(1)**2+DKxyz(2)**2))*(SIN(DUMMY2)/DUMMY2)**2).LT.DOWN_INT) THEN THETA_SIGNAL_STOP1=TEMP ELSE THETA_SIGNAL_START1=TEMP ENDIF ENDDO THETA_STOP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2 THETA_MAX=MAX(THETA_MAX,THETA_STOP) THETA_SIGNAL_START1=0 THETA_SIGNAL_STOP1=ANG1N DO WHILE ((abs(THETA_SIGNAL_STOP1-THETA_SIGNAL_START1)).GT.TOLERANCE) TEMP=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2. THETA(1) = TEMP CALL DeltaK(N,THETA,DELTA_K) DUMMY2=0.5*CryLength*DKxyz(3) IF ((EXP(-0.5*Pump_W**2*(DKxyz(1)**2+DKxyz(2)**2))*(SIN(DUMMY2)/DUMMY2)**2).LT.DOWN_INT) THEN THETA_SIGNAL_START1=TEMP ELSE THETA_SIGNAL_STOP1=TEMP ENDIF ENDDO THETA_START=(THETA_SIGNAL_START1+THETA_SIGNAL_STOP1)/2 THETA_MIN=MIN(THETA_MIN,THETA_START) ITER1=ITER1-INC ANG1N=(THETA_START+THETA_STOP)/2 ENDDO ! Conversion Internal angle -> External angle ! ANG1=abs(ANG1) ANG1ext = asin(sin(ANG1)*N_SIG) DT_SIG_P_ext = asin(sin(THETA_MAX)*N_SIG) DT_SIG_N_ext = asin(sin(THETA_MIN)*N_SIG) WRITE(OUTCHAN,"(F10.5,',',F20.7,',',F10.7,',',F10.3)") LAMBDA_SIGNAL,THETA_MAX*180/PI,DT_SIG_P_ext*180/PI,DOWN_INT WRITE(OUTCHAN,"(F10.5,',',F20.7,',',F10.7,',',F10.3)") LAMBDA_SIGNAL,THETA_MIN*180/PI,DT_SIG_N_ext*180/PI,DOWN_INT ENDDO ENDIF ENDIF ENDDO CLOSE(INCHAN) CLOSE(OUTCHAN) pause !*************************************************************************** ELSE IF ( MENU .EQ. 10 ) THEN ! CALL OutwindowClear WRITE(*, *) '(10) 3D plot, Phase-matching function = f(Phi,Lambda)' WRITE(OUTCHAN, *) '(10) 3D plot, Phase-matching function = f(Phi,Lambda)' WRITE(*, "(' Wavelength of Pump (microns) =', $)") READ(INCHAN, *) LAMBDA_PUMP WRITE(OUTCHAN,"(' Wavelength of Pump (microns) =', F10.5)") LAMBDA_PUMP WRITE(OUTCHAN,*) ' Nx, Ny, Nz' CALL INDEX(CRYSTAL, LAMBDA_PUMP, NX,NY, NZ, ABS_MIN, ABS_MAX) WRITE(OUTCHAN,"('Pump: ',F10.5,F10.5,F10.5)") NX,NY,NZ WRITE(*,"(' Phi Pump (degrees) =',$)") READ(INCHAN,*) ALPHA WRITE(OUTCHAN,"(' Phi Pump (deg)=',F10.5)") ALPHA WRITE(*,"(' Theta Pump (degrees) =',$)") READ(INCHAN,*) BETA WRITE(OUTCHAN,"(' Theta Pump (deg)=',F10.5)") BETA WRITE(*,"(' Phi Signal (degrees) =',$)") READ(INCHAN,*) PHI_SIG2 WRITE(OUTCHAN,"(' Phi Signal (deg)=',F10.5)") PHI_SIG2 ALPHA = ALPHA * PI/ 180. BETA = BETA * PI/ 180. PHI_SIG2= PHI_SIG2 * PI /180. WRITE(*,*) 'Phase Matching TYPEPMs' WRITE(*,*) ' TYPEPM Pump Signal Idler' WRITE(*,*) ' 1 FAST SLOW SLOW' WRITE(*,*) ' 2 FAST FAST SLOW' WRITE(*,*) ' 3 FAST SLOW FAST' WRITE(*,"(' Choice: ',$)") READ(INCHAN,*) TYPEPM IF ( TYPEPM .EQ. 1 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Slow Idler=Slow' ELSE IF ( TYPEPM .EQ. 2 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Fast Idler=Slow' ELSE IF ( TYPEPM .EQ. 3 ) THEN WRITE(OUTCHAN,*) 'Pump=Fast Signal=Slow Idler=Fast' ENDIF WRITE(*,"(' Length of Crystal (microns) =',$)") READ(INCHAN,*) CryLength WRITE(OUTCHAN,"(' Crystal length =',F10.0)") CryLength WRITE(*, "(' Gaussian Pump (intensity FWHM in microns) =', $)") READ(INCHAN, *) Pump_FWHM WRITE(OUTCHAN, "(' Gaussian Pump FWHM =', F10.0)") Pump_FWHM Pump_W=Pump_FWHM*0.8493218 WRITE(*, "(' Lambda Start (microns) =', $)") READ(INCHAN, *) LAMBDA_START WRITE(OUTCHAN, "(' Lambda Start (microns) =', F10.5)") LAMBDA_START WRITE(*, "(' Lambda Stop (microns) =', $)") READ(INCHAN, *) LAMBDA_STOP WRITE(OUTCHAN, "(' Lambda Stop (microns) =', F10.5)") LAMBDA_STOP WRITE(*, "(' Downconversion Intensity inc =', $)") READ(INCHAN, *) DOWN_INT_INC WRITE(OUTCHAN, "(' Downconversion Intensity inc =', F10.3)") DOWN_INT_INC CALL KPUMP PHI_SIGNAL=PHI_SIG2 PHI_IDLER=PHI_SIGNAL+PI WRITE(OUTCHAN,*) 'Lambda Signal',',','Phi Signal',',','Downconversion Intensity' ANG1=XINIT1*FACTOR ANG2=XINIT2*FACTOR WRITE(*,*) 'Lambda signal = ','Theta signal optimum = ' F=1 DO LAMBDA_SIGNAL =LAMBDA_START,LAMBDA_STOP,(LAMBDA_STOP-LAMBDA_START)/100 LAMBDA_IDLER = 1/(1/LAMBDA_PUMP - 1/LAMBDA_SIGNAL) IF ( LAMBDA_SIGNAL .GE. ABS_MIN .AND. & LAMBDA_SIGNAL .LE. ABS_MAX .AND. & LAMBDA_PUMP .GE. ABS_MIN .AND. & LAMBDA_PUMP .LE. ABS_MAX .AND. & LAMBDA_IDLER .GE. ABS_MIN .AND. & LAMBDA_IDLER .LE. ABS_MAX ) THEN PHI_SIGNAL=PHI_SIG2 PHI_IDLER=PHI_SIGNAL+PI LW=N*(N+10) ! Initial guesses for theta (signal,idler) optimum IF ((MENU3==2).AND.(F<0.00001)) THEN X0(1)=ANG1 X0(2)=ANG2 ELSE X0(1)=XINIT1*FACTOR X0(2)=XINIT2*FACTOR ENDIF CALL UNCMND (N,X0,DeltaK,THETA,F,INFO,W,LW) ! THETA(1)=ABS(THETA(1)) ANG1=THETA(1) ! THETA(2)=ABS(THETA(2)) ANG2=THETA( 2 ) IF (F<0.00001) THEN WRITE(*,"(F10.5,F10.5)") LAMBDA_SIGNAL,THETA(1)*180/PI DO DOWN_INT= 1. , 0.1 , -DOWN_INT_INC PHI_SIGNAL=PHI_SIG2 PHI_IDLER=PHI_SIG2+PI PHI_START=PHI_SIG2 PHI_STOP=PHI_SIG2+PI/2 DO WHILE ((PHI_STOP-PHI_START).GT.TOLERANCE) TEMP=(PHI_START+PHI_STOP)/2. PHI_SIGNAL = TEMP CALL DeltaK(N,THETA,DELTA_K) DUMMY2=0.5*CryLength*DKxyz(3) IF ((EXP(-0.5*Pump_W**2*(DKxyz(1)**2+DKxyz(2)**2))*(SIN(DUMMY2)/DUMMY2)**2).LT.DOWN_INT) THEN PHI_STOP=TEMP ELSE PHI_START=TEMP ENDIF ENDDO PHI_MAX=(PHI_START+PHI_STOP)/2 PHI_MIN=2*PHI_SIG2-PHI_MAX WRITE(OUTCHAN,"(F10.5,',',F20.7,',',F10.3)") LAMBDA_SIGNAL,PHI_MAX*180/PI,DOWN_INT WRITE(OUTCHAN,"(F10.5,',',F20.7,',',F10.3)") LAMBDA_SIGNAL,PHI_MIN*180/PI,DOWN_INT ENDDO ENDIF ENDIF ENDDO CLOSE(INCHAN) CLOSE(OUTCHAN) pause !*************************************************************************** ELSE IF ( MENU .EQ. 11) THEN ! CALL OutwindowClear WRITE(*, *) "" WRITE(*, *) 'HELP HELP HELP HELP HELP' WRITE(*, *) "" WRITE(*, *) 'The main menu gives a list of 11 options. Within each of the options, other' WRITE(*, *) 'information such as a file name, crystal length or TYPEPM, wavelengths or' WRITE(*, *) 'angles of interest and other information may be needed. After an operation' WRITE(*, *) 'is complete, the program will return to the main menu and another option ' WRITE(*, *) 'may be entered.' WRITE(*, *) "" WRITE(*, *) 'This program does not plot data. It creates output text files in a table' WRITE(*, *) 'format, with entries delimited by commas (,) to be imported into suitable' WRITE(*, *) 'plotting applications.' WRITE(*, *) "" WRITE(*, *) 'For further help, contact Alan Migdall or David Branning' WRITE(*, *) '' WRITE(*, *) "" WRITE(*, *) 'The available options are:' pause ! CALL OutwindowClear WRITE(*, *) '(1) 2D Plot, (Nx,Ny,Nz) = f(Lambda)' WRITE(*, *) 'This selection allows the indices of refraction for the three directions' WRITE(*, *) 'to be viewed as a function of wavelength, giving insight into which' WRITE(*, *) 'wavelengths are suitable for phase matching. The output culumns are four columns' WRITE(*, *) '(three in the case of uniaxial crystals) containing the indices of' WRITE(*, *) 'refraction versus wavelength.' WRITE(*, *) "" WRITE(*, *) 'The input may come from a text file named PHASEMATCHIN1 formatted as follows:' WRITE(*, *) "" WRITE(*, *) 'Output file name' WRITE(*, *) 'Crytal' WRITE(*, *) "" WRITE(*, *) 'The output columns are : Lambda, Nx,Ny,Nz' pause ! CALL OutwindowClear WRITE(*, *) '(2) 3D Plot, nSlow-nFast = f(Theta,Phi)' WRITE(*, *) 'This option is intended to help in designing a suitable orientation for' WRITE(*, *) 'the crystal. The difference in indices of refraction, n-Slow - n-Fast,' WRITE(*, *) 'is given as output in terms of the direction of the pump through the ' WRITE(*, *) 'crystal. Two angles, Theta and Phi, specify this direction relative to' WRITE(*, *) 'the optic axis (only theta is necessary for uniaxial crystals). The ' WRITE(*, *) 'output should be plotted in three dimensions (the output could be displayed' WRITE(*, *) 'by shading or color rather than a third dimension). The difference in' WRITE(*, *) 'slow and fast indices is given at the pump wavelength; large differences in' WRITE(*, *) '(n-Slow - n-Fast) correspond to regions where phase-matching is probable.' WRITE(*, *) "" WRITE(*, *) 'The input may come from a text file named PHASEMATCHIN2 formatted as follows:' WRITE(*, *) "" WRITE(*, *) 'Output file name' WRITE(*, *) 'Crytal' WRITE(*, *) 'Scale factor' WRITE(*, *) 'Lambda pump' WRITE(*, *) 'Phi pump inc' WRITE(*, *) 'Theta pump inc' WRITE(*, *) "" WRITE(*, *) 'The output culumns are : Theta, Phi, nSlow-nFast' pause ! CALL OutwindowClear WRITE(*, *) '(3) 3D Plot, minimum Delta-K = f(Theta,Phi)' WRITE(*, *) 'Delta-K is the difference between the pump propagation vector and the sum' WRITE(*, *) 'of the signal and idler propagation vectors:' WRITE(*, *) "" WRITE(*, *) 'Delta-K = Kpump - (Ksignal + Kidler), with LambdaPump and LambdaSignal given' WRITE(*, *) "" WRITE(*, *) 'As an alternative to option 3, this selection may help in determining an' WRITE(*, *) 'appropriate crystal orientation. The output is given in the same format as' WRITE(*, *) 'in option 3, except that the minimum value of Delta-K is given instead of' WRITE(*, *) 'the quantity (n-Slow - n-Fast). Phase matching can occur in infinitely long' WRITE(*, *) 'crystals only for the orientations that result in Delta-K = 0. For crystals' WRITE(*, *) 'of finite length, downconversion may still occur for small nonzero values ' WRITE(*, *) 'of Delta-K. ' WRITE(*, *) "" WRITE(*, *) 'The input may come from a text file named PHASEMATCHIN3 formatted as follows:' WRITE(*, *) "" WRITE(*, *) 'Output file name' WRITE(*, *) 'Crytal' WRITE(*, *) 'Scale factor' WRITE(*, *) 'Lambda pump' WRITE(*, *) 'Lambda signal' WRITE(*, *) 'Phi pump inc' WRITE(*, *) 'Theta pump inc' WRITE(*, *) 'TYPEPM' WRITE(*, *) "" WRITE(*, *) 'The output culumns are : Theta, Phi, Delta-K' pause ! CALL OutwindowClear WRITE(*, *) '(4) 3D Plot, Phase-matching function = P(Dk-transverse,Dkz) ' WRITE(*, *) 'The downconversion process is governed by the phase-matching function (P),' WRITE(*, *) 'given by:' WRITE(*, *) '' WRITE(*, *) 'P = exp(-.5*W^2*Dktransverse^2)*(sin(.5*L*Dkz)/(.5*L*Dkz)^2' WRITE(*, *) '(Crystal length and pump width given) ' WRITE(*, *) '' WRITE(*, *) 'When P is one, downconversion will be at a maximum; no downconversion may' WRITE(*, *) 'occur for pump, signal, and idler combinations yielding P = 0. For crystals ' WRITE(*, *) 'of finite length and pump beams of finite width, P does not fall ' WRITE(*, *) 'immediately from one to zero for nonzero values of Delta-K. This option ' WRITE(*, *) 'gives the values of P as Delta-K is varied along the pump beam (in the' WRITE(*, *) 'z direction) and transverse to this direction, showing the region of ' WRITE(*, *) 'Delta-K mismatch which still allows some downconversion to occur. The ' WRITE(*, *) 'output can be a 3D plot or a collection of curves which represent ' WRITE(*, *) 'different fixed values of the Phase-matching function. ' WRITE(*, *) "" WRITE(*, *) 'The input may come from a text file named PHASEMATCHIN4 formatted as follows:' WRITE(*, *) "" WRITE(*, *) 'Crytal length' WRITE(*, *) 'Pump width' WRITE(*, *) "" WRITE(*, *) 'For a range of Ps from 1 to 0.1 the output culumns are: DK-transverse,DKz,P' pause ! CALL OutwindowClear WRITE(*, *) '(5) 2D Plot, Theta signal vs. Theta idler (lambda pump and lambda signal fixed)' WRITE(*, *) 'This option is related to option 5, but calculates the range of signal and ' WRITE(*, *) 'idler angles which result in a decrease of the phase-matching function (P)' WRITE(*, *) 'from 1 to some user-defined value. For example, if P = .5 is chosen,' WRITE(*, *) 'the reported range of angles for theta signal will correspond to the' WRITE(*, *) 'observed FWHM of the downconversion.' WRITE(*, *) "" WRITE(*, *) 'The input may come from a text file named PHASEMATCHIN5 formatted as follows:' WRITE(*, *) "" WRITE(*, *) 'Output file name' WRITE(*, *) 'Crytal' WRITE(*, *) 'Scale factor' WRITE(*, *) 'Lambda pump' WRITE(*, *) 'Lambda signal' WRITE(*, *) 'Phi pump' WRITE(*, *) 'Theta pump' WRITE(*, *) 'Phi signal' WRITE(*, *) 'TYPEPM' WRITE(*, *) 'Crytal length' WRITE(*, *) 'Pump width' WRITE(*, *) 'Downconversion intensity' WRITE(*, *) 'Theta inc' WRITE(*, *) "" WRITE(*, *) 'The output culumns are (ThetaSignal, ThetaIdler) points mapping out where' WRITE(*, *) 'P=the selected value' pause ! CALL OutwindowClear WRITE(*, *) '(6) Polar plot, (Optimum Theta Idler,Optimum Theta Signal) = f(Phi)' WRITE(*, *) 'Once a crystal orientation is chosen, this option allows the output at a' WRITE(*, *) 'specific wavelength to be examined. Given a fixed signal output wavelength,' WRITE(*, *) 'external angles for both it and the correlated idler wavelength are' WRITE(*, *) 'returned. The external angles are given in polar coordinates with theta' WRITE(*, *) 'corresponding to the cone opening angle from the pump direction and phi' WRITE(*, *) 'corresponding to the azimuthal angle about the pump direction :' WRITE(*, *) '"phi=0 degrees" being in the plane formed by the pump direction and crystal' WRITE(*, *) 'axis (uniaxial case) or in the plane containing the crystal axes' WRITE(*, *) '(biaxial case)' WRITE(*, *) "" WRITE(*, *) 'The input may come from a text file named PHASEMATCHIN6 formatted as follows:' WRITE(*, *) "" WRITE(*, *) 'Output file name' WRITE(*, *) 'Crytal' WRITE(*, *) 'Scale factor' WRITE(*, *) 'Lambda pump' WRITE(*, *) 'Lambda signal' WRITE(*, *) 'Phi pump' WRITE(*, *) 'Theta pump' WRITE(*, *) 'Phi signal inc' WRITE(*, *) 'TYPEPM' WRITE(*, *) "" WRITE(*, *) 'The output culumns are : phi, theta signal internal, theta idler internal,' WRITE(*, *) ' theta signal external, theta idler external' pause ! CALL OutwindowClear WRITE(*, *) '(7) and (8) 2D Plot, Optimum Theta = f(lambda) at chosen Phi' WRITE(*, *) ' with spreading in Theta,Phi' WRITE(*, *) 'Once a crystal, pump wavelength, and phase-matching TYPEPM have been chosen,' WRITE(*, *) 'this selection gives the optimum angle (theta signal) that results in perfect' WRITE(*, *) 'phase-matching (P = 1 or delta-K = 0) as a function of wavelength. No output is ' WRITE(*, *) 'generated for wavelengths that cannot produce perfect phase-matching.' WRITE(*, *) "" WRITE(*, *) 'Also, because of the finite crystal length and non-zero pump width,' WRITE(*, *) 'imperfect phase-matching can still occur if theta is not equal to theta-' WRITE(*, *) 'optimum and if phi is not equal to the chosen phi. Therefore,this option' WRITE(*, *) 'also reports the spread in the allowed values of theta and phi' WRITE(*, *) 'which result in a decrease of the phase-matching function to some value' WRITE(*, *) 'chosen by the user. For example, if P = .5 is chosen, the reported spreads' WRITE(*, *) 'in theta and phi will correspond to the FWHM of the downconversion intensity.' WRITE(*, *) 'In option 7 Theta idler is fixed to Theta idler optimum.' WRITE(*, *) 'In option 8 Theta idler is allowed to change.' WRITE(*, *) "" WRITE(*, *) 'The input may come from a text file named PHASEMATCHIN7(8) formatted as' WRITE(*, *) 'follows:' WRITE(*, *) "" WRITE(*, *) 'Output file name' WRITE(*, *) 'Crytal' WRITE(*, *) 'Scale factor' WRITE(*, *) 'Lambda pump' WRITE(*, *) 'Phi pump' WRITE(*, *) 'Theta pump' WRITE(*, *) 'Phi signal' WRITE(*, *) 'TYPEPM' WRITE(*, *) 'Crytal length' WRITE(*, *) 'Pump width' WRITE(*, *) 'Lambda start' WRITE(*, *) 'Lambda stop' WRITE(*, *) 'Downconversion intensity' WRITE(*, *) 'Theta inc (Option 8 only)' WRITE(*, *) "" WRITE(*, *) 'The output culumns are : Lambda,ThetaSignal_Optimum Internal,' WRITE(*, *) ' ThetaSignal Spread+Internal,Theta Signal Spread-Internal,' WRITE(*, *) ' ThetaSignal_Optimum External,ThetaSignal Spread+External,' WRITE(*, *) ' Theta Signal Spread-External,' WRITE(*, *) ' Phi Spread+,Phi Spread-' pause ! CALL OutwindowClear WRITE(*, *) '(9) and (10) 3D Plot, Phase-matching function = P(Theta,Phi,Lambda)' WRITE(*, *) 'These options give the same result as option 8, but also computes the' WRITE(*, *) 'in theta and phi corresponding to many different values of the spreads ' WRITE(*, *) 'Phase-matching function between 0.1 and 1. A three-dimensional plot of' WRITE(*, *) 'this data will generate a picture of the emitted downconversion intensity' WRITE(*, *) 'as a function of wavelength and theta (9) or wavelength and phi (10).' WRITE(*, *) 'Again, no output is generated for wavelengths that cannot produce perfect' WRITE(*, *) 'phase-matching (P = 1).' WRITE(*, *) "" WRITE(*, *) 'The input may come from a text file named PHASEMATCHIN9(10) formatted as' WRITE(*, *) 'follows:' WRITE(*, *) "" WRITE(*, *) 'Output file name' WRITE(*, *) 'Crytal' WRITE(*, *) 'Scale factor' WRITE(*, *) 'Lambda pump' WRITE(*, *) 'Phi pump' WRITE(*, *) 'Theta pump' WRITE(*, *) 'Phi signal' WRITE(*, *) 'TYPEPM' WRITE(*, *) 'Crytal length' WRITE(*, *) 'Pump width' WRITE(*, *) 'Lambda start' WRITE(*, *) 'Lambda stop' WRITE(*, *) 'Downconversion intensity inc' WRITE(*, *) 'Theta inc (Option 9 only)' WRITE(*, *) "" WRITE(*, *) 'The output culumns are : Lambda,Theta,Downconversion intensity (9)' WRITE(*, *) 'The output culumns are : Lambda,Phi,Downconversion intensity (10)' pause ! CALL OutwindowClear WRITE(*, *) 'HINTS for running options 7-10.' WRITE(*, *) 'This program uses a iterative minimization routine that is sensitive to' WRITE(*, *) 'an initial guess for downconversion opening angle THETA.' WRITE(*, *) 'An initial guess is provided, but may be changed by a multiplier.' WRITE(*, *) "" WRITE(*, *) 'An additional option that may be chosen is whether to use the same initial' WRITE(*, *) 'guess for the minimization routine at each wavelength, or, after the solution' WRITE(*, *) 'or the first wavelength, to use the latest result as the initial guess for' WRITE(*, *) 'the opening angle of the next wavelength. The most appropriate option will' WRITE(*, *) 'depend on the specific problem being solved. The first option is useful to' WRITE(*, *) 'reveal multiple THETA solutions for a given PHI, while the second option will' WRITE(*, *) 'tend to follow whichever solution is found initially.' pause ! CALL OutwindowClear WRITE(*, *) '(12) Stop' WRITE(*, *) 'This exits the program.' pause ENDIF ENDDO END !*************************************************************************** !****** FUNCTIONS AND SUBROUTINES ************************************ ! This subroutine calculates delta k for parameters passed mostly by common block SUBROUTINE DeltaK(N,THETA,F) IMPLICIT NONE REAL*8 F,PI ! TYPEPM of crystal for phase matching INTEGER CRYSTAL,N ! TYPEPM of phase matching INTEGER TYPEPM ! Pump phi direction in crystal coordinate system REAL*8 ALPHA ! Pump theta direction in crystal coordinate system REAL*8 BETA REAL*8 PP2,QQ4 ! wavelength of idler,pump,signal REAL*8 LAMBDA_IDLER,LAMBDA_PUMP,LAMBDA_SIGNAL REAL*8 ABS_MIN,ABS_MAX REAL*8 N_IDLER,N_PUMP,N_SIG REAL*8 NX_IDLER,NX_SIG,NY_IDLER,NY_SIG,NZ_IDLER,NZ_SIG ! phi of idler in pump coordinate system REAL*8 PHI_IDLER ! phi of signal in pump coordinate system REAL*8 PHI_SIGNAL REAL*8 DKxyz(3) REAL*8 PP_IDLER,PP_SIG,QQ_IDLER,QQ_SIG ! theta of signal (1) and idler (2) in pump coordinate system REAL*8 THETA(2) REAL*8 TX_IDLER,TX_PUMP,TX_SIG REAL*8 TY_IDLER,TY_PUMP,TY_SIG REAL*8 TZ_IDLER,TZ_PUMP,TZ_SIG PARAMETER (PI = 3.141592653589793238) COMMON /CM/ CRYSTAL, TYPEPM, TX_PUMP,TY_PUMP,TZ_PUMP, & N_PUMP, N_SIG, N_IDLER, ALPHA, BETA, PHI_IDLER, PHI_SIGNAL, & LAMBDA_IDLER,LAMBDA_SIGNAL,LAMBDA_PUMP,DKxyz ! TX,TY,TZ are in crystal coordinates IF ( TYPEPM .EQ. 1 .OR. TYPEPM .EQ. 2 .OR. TYPEPM .EQ. 3 ) THEN CALL INDEX(CRYSTAL, LAMBDA_SIGNAL, NX_SIG,NY_SIG, NZ_SIG, ABS_MIN, ABS_MAX) TX_SIG = COS(ALPHA)*SIN(BETA)*COS(THETA(1)) + & COS(ALPHA)*COS(BETA)*SIN(THETA(1))*COS(PHI_SIGNAL) - & SIN(ALPHA)*SIN(THETA(1))*SIN(PHI_SIGNAL) TY_SIG = SIN(ALPHA)*SIN(BETA)*COS(THETA(1)) + & SIN(ALPHA)*COS(BETA)*SIN(THETA(1))*COS(PHI_SIGNAL) + & COS(ALPHA)*SIN(THETA(1))*SIN(PHI_SIGNAL) TZ_SIG = COS(BETA)*COS(THETA(1)) - & SIN(BETA)*SIN(THETA(1))*COS(PHI_SIGNAL) ! The index of refraction PP_SIG = (TX_SIG**2)*(1/NY_SIG**2 + 1/NZ_SIG**2) + & (TY_SIG**2)*(1/NX_SIG**2 + 1/NZ_SIG**2) + & (TZ_SIG**2)*(1/NX_SIG**2 + 1/NY_SIG**2) QQ_SIG = (TX_SIG**2)/(NY_SIG**2 * NZ_SIG**2) + & (TY_SIG**2)/(NX_SIG**2 * NZ_SIG**2) + & (TZ_SIG**2)/(NX_SIG**2 * NY_SIG**2) CALL INDEX(CRYSTAL, LAMBDA_IDLER, NX_IDLER, & NY_IDLER, NZ_IDLER, ABS_MIN, ABS_MAX) ! Direction cosines, with respect to crytal axes system, with phi ! and theta in pump axes system. TX_IDLER = COS(ALPHA)*SIN(BETA)*COS(THETA(2)) + & COS(ALPHA)*COS(BETA)*SIN(THETA(2))*COS(PHI_IDLER) - & SIN(ALPHA)*SIN(THETA(2))*SIN(PHI_IDLER) TY_IDLER = SIN(ALPHA)*SIN(BETA)*COS(THETA(2)) + & SIN(ALPHA)*COS(BETA)*SIN(THETA(2))*COS(PHI_IDLER) + & COS(ALPHA)*SIN(THETA(2))*SIN(PHI_IDLER) TZ_IDLER = COS(BETA)*COS(THETA(2)) - & SIN(BETA)*SIN(THETA(2))*COS(PHI_IDLER) ! The index of refraction PP_IDLER = (TX_IDLER**2)*(1/NY_IDLER**2 + 1/NZ_IDLER**2) + & (TY_IDLER**2)*(1/NX_IDLER**2 + 1/NZ_IDLER**2) + & (TZ_IDLER**2)*(1/NX_IDLER**2 + 1/NY_IDLER**2) QQ_IDLER = (TX_IDLER**2)/(NY_IDLER**2 * NZ_IDLER**2) + & (TY_IDLER**2)/(NX_IDLER**2 * NZ_IDLER**2) + & (TZ_IDLER**2)/(NX_IDLER**2 * NY_IDLER**2) ! Now get the fast or slow. IF ( TYPEPM .EQ. 1 ) THEN N_SIG = SQRT( 2.0/(PP_SIG - SQRT(PP_SIG**2 - 4*QQ_SIG)) ) N_IDLER = SQRT( 2.0/(PP_IDLER- & SQRT(PP_IDLER**2 - 4*QQ_IDLER)) ) ELSE IF ( TYPEPM .EQ. 2 ) THEN N_SIG = SQRT( 2.0/(PP_SIG + SQRT(PP_SIG**2 - 4*QQ_SIG)) ) N_IDLER = SQRT( 2.0/(PP_IDLER - & SQRT(PP_IDLER**2 - 4*QQ_IDLER)) ) ELSE IF ( TYPEPM .EQ. 3 ) THEN N_SIG = SQRT( 2.0/(PP_SIG - SQRT(PP_SIG**2 - 4*QQ_SIG)) ) N_IDLER = SQRT( 2.0/(PP_IDLER + & SQRT(PP_IDLER**2 - 4*QQ_IDLER)) ) ENDIF ! Now convert TX, TY, TZ back to pump (") coordinates TX_SIG=SIN(THETA(1))*COS(PHI_SIGNAL) TY_SIG=SIN(THETA(1))*SIN(PHI_SIGNAL) TZ_SIG=COS(THETA(1)) TX_IDLER=SIN(THETA(2))*COS(PHI_IDLER) TY_IDLER=SIN(THETA(2))*SIN(PHI_IDLER) TZ_IDLER=COS(THETA(2)) DKxyz(1) = 2*PI *(N_SIG*TX_SIG/LAMBDA_SIGNAL + N_IDLER*TX_IDLER/LAMBDA_IDLER) DKxyz(2) = 2*PI *(N_SIG*TY_SIG/LAMBDA_SIGNAL + N_IDLER*TY_IDLER/LAMBDA_IDLER) DKxyz(3) = 2*PI *(N_SIG*TZ_SIG/LAMBDA_SIGNAL + N_IDLER*TZ_IDLER/LAMBDA_IDLER & - N_PUMP/LAMBDA_PUMP) F= DKxyz(1)**2+DKxyz(2)**2+DKxyz(3)**2 F = SQRT(F) ELSE CALL ERR_MSG("No such TYPEPM", 13) ! End of crystal TYPEPM. ENDIF ! RETURN END !*************************************************************************** ! Subroutine: KPUMP ! Initializes N_PUMP TX_PUMP TY_PUMP and TZ_PUMP ! The index is the Fast index SUBROUTINE KPUMP INTEGER CRYSTAL,TYPEPM REAL*8 NX,NY,NZ REAL*8 ABS_MIN, ABS_MAX REAL*8 TX_PUMP,TY_PUMP,TZ_PUMP REAL*8 N_PUMP,N_SIG,N_IDLER REAL*8 ALPHA,BETA REAL*8 PS,QQ REAL*8 PHI_IDLER,PHI_SIGNAL REAL*8 LAMBDA_IDLER,LAMBDA_SIGNAL,LAMBDA_PUMP,DKxyz(3) COMMON /CM/ CRYSTAL, TYPEPM, TX_PUMP,TY_PUMP,TZ_PUMP, & N_PUMP, N_SIG, N_IDLER, ALPHA, BETA, PHI_IDLER, PHI_SIGNAL, & LAMBDA_IDLER,LAMBDA_SIGNAL,LAMBDA_PUMP , DKxyz TX_PUMP = COS(ALPHA)*SIN(BETA) TY_PUMP = SIN(ALPHA)*SIN(BETA) TZ_PUMP = COS(BETA) CALL INDEX(CRYSTAL, LAMBDA_PUMP, NX,NY, NZ, ABS_MIN, ABS_MAX) PS = (TX_PUMP**2)*(1/NY**2 + 1/NZ**2) + & (TY_PUMP**2)*(1/NX**2 + 1/NZ**2) + & (TZ_PUMP**2)*(1/NX**2 + 1/NY**2) QQ = (TX_PUMP**2)/(NY**2 * NZ**2) + & (TY_PUMP**2)/(NX**2 * NZ**2) + & (TZ_PUMP**2)/(NX**2 * NY**2) N_PUMP = SQRT( 2.0/(PS + SQRT(PS**2 - 4*QQ)) ) END !*************************************************************************** SUBROUTINE ERR_MSG(MESSAGE, LEN) ! SUBPROGRAM: ERR_MSG ! AUTHOR: Jeffrey S. Orszak ! DATE CREATED: October 1995 ! MODIFIED: 22 February 1996 ! COPYRIGHT: This work is written under contract to the ! National Institute of Standards and Technology, as such ! it is not subject to U. S. copyright. ! Routine for error message handling. IMPLICIT NONE ! Error message CHARACTER*132 MESSAGE ! Internal format statement CHARACTER*7 FNAME ! Length of error message. INTEGER LEN WRITE(FNAME, 100) LEN 100 FORMAT( "(A", I2, ")") ! Should actually write to STDERR WRITE (*, FNAME) MESSAGE STOP "Error" END !*************************************************************************** SUBROUTINE LIST IMPLICIT NONE WRITE(*, *) 'Uniaxial crystals ' WRITE(*, *) '(Nx = No, Ny = No, Nz = Ne)' WRITE(*, *) ' 10 Ag3AsS3 ( proustite )' WRITE(*, *) ' 20 AgGaSe2 Ref 1' WRITE(*, *) ' 30 AgGaSe2 Ref 2' WRITE(*, *) ' 40 BBO Ref 1' WRITE(*, *) ' 50 BBO Ref 2' WRITE(*, *) ' 60 GaSe' WRITE(*, *) ' 70 KDP Ref 1' WRITE(*, *) ' 80 KDP Ref 2' WRITE(*, *) ' 90 KD*P Ref 1' WRITE(*, *) ' 100 KD*P Ref 2' WRITE(*, *) ' 110 K3Li2Nb5O15' WRITE(*, *) ' 120 LiIO3 Ref 6 (use Ref. 7)' WRITE(*, *) ' 130 LiIO3 Ref 7' WRITE(*, *) ' 140 LiNO3' WRITE(*, *) ' 150 Tl3AsSe3' WRITE(*, *) ' 160 ZnGeP2' WRITE(*, *) ' 170 AgGaS2' WRITE(*, *) 'Biaxial crystals ' WRITE(*, *) ' 200 KTP' WRITE(*, *) ' 210 KTA' WRITE(*, *) ' 220 Magnesium-Barium Flouride' WRITE(*, *) ' 230 KNbO3 ( 22 C )' WRITE(*, *) ' 231 KNbO3 ( 50 C )' WRITE(*, *) ' 232 KNbO3 ( 75 C )' WRITE(*, *) ' 233 KNbO3 (100 C )' WRITE(*, *) ' 234 KNbO3 (140 C )' WRITE(*, *) ' 235 KNbO3 (180 C )' END !*************************************************************************** SUBROUTINE FLIST(OUTCHAN,CRYSTAL) IMPLICIT NONE INTEGER OUTCHAN INTEGER CRYSTAL REAL*8 LAMBDA REAL*8 ABS_MIN,ABS_MAX REAL*8 NX,NY,NZ IF ( CRYSTAL .EQ. 10) THEN WRITE(OUTCHAN, *) 'Crystal: Ag3AsS3 ( proustite )' ELSE IF ( CRYSTAL .EQ. 20) THEN WRITE(OUTCHAN, *) 'Crystal: AgGaSe2 Ref 1' ELSE IF ( CRYSTAL .EQ. 30) THEN WRITE(OUTCHAN, *) 'Crystal: AgGaSe2 Ref 2' ELSE IF ( CRYSTAL .EQ. 40) THEN WRITE(OUTCHAN, *) 'Crystal: BBO Ref 1' ELSE IF ( CRYSTAL .EQ. 50) THEN WRITE(OUTCHAN, *) 'Crystal: BBO Ref 2' ELSE IF ( CRYSTAL .EQ. 60) THEN WRITE(OUTCHAN, *) 'Crystal: GaSe' ELSE IF ( CRYSTAL .EQ. 70) THEN WRITE(OUTCHAN, *) 'Crystal: KDP Ref 1' ELSE IF ( CRYSTAL .EQ. 80) THEN WRITE(OUTCHAN, *) 'Crystal: KDP Ref 2' ELSE IF ( CRYSTAL .EQ. 90) THEN WRITE(OUTCHAN, *) 'Crystal: KD*P Ref 1' ELSE IF ( CRYSTAL .EQ. 100) THEN WRITE(OUTCHAN, *) 'Crystal: KD*P Ref 2' ELSE IF ( CRYSTAL .EQ. 110) THEN WRITE(OUTCHAN, *) 'Crystal: K3Li2Nb5O15' ELSE IF ( CRYSTAL .EQ. 120) THEN WRITE(OUTCHAN, *) 'Crystal: LiIO3 Ref 6 (use Ref. 7)' ELSE IF ( CRYSTAL .EQ. 130) THEN WRITE(OUTCHAN, *) 'Crystal: LiIO3 Ref 7' ELSE IF ( CRYSTAL .EQ. 140) THEN WRITE(OUTCHAN, *) 'Crystal: LiNO3' ELSE IF ( CRYSTAL .EQ. 150) THEN WRITE(OUTCHAN, *) 'Crystal: Tl3AsSe3' ELSE IF ( CRYSTAL .EQ. 160) THEN WRITE(OUTCHAN, *) 'Crystal: ZnGeP2' ELSE IF ( CRYSTAL .EQ. 170) THEN WRITE(OUTCHAN, *) 'Crystal: AgGaS2' ELSE IF ( CRYSTAL .EQ. 200) THEN WRITE(OUTCHAN, *) 'Crystal: KTP' ELSE IF ( CRYSTAL .EQ. 210) THEN WRITE(OUTCHAN, *) 'Crystal: KTA' ELSE IF ( CRYSTAL .EQ. 220) THEN WRITE(OUTCHAN, *) 'Crystal: Magnesium-Barium Flouride' ELSE IF ( CRYSTAL .EQ. 230) THEN WRITE(OUTCHAN, *) 'Crystal: KNbO3 ( 22 C )' WRITE(OUTCHAN, *) 'J. Opt. Soc. Am. B/Vol.9,No.3/March 1992' WRITE(OUTCHAN, *) 'Refractive indices of orthorhombic KNbO3' WRITE(OUTCHAN, *) 'B. Zysset' ELSE IF ( CRYSTAL .EQ. 231) THEN WRITE(OUTCHAN, *) 'Crystal: KNbO3 ( 50 C )' WRITE(OUTCHAN, *) 'J. Opt. Soc. Am. B/Vol.9,No.3/March 1992' WRITE(OUTCHAN, *) 'Refractive indices of orthorhombic KNbO3' WRITE(OUTCHAN, *) 'B. Zysset' ELSE IF ( CRYSTAL .EQ. 232) THEN WRITE(OUTCHAN, *) 'Crystal: KNbO3 ( 75 C )' WRITE(OUTCHAN, *) 'J. Opt. Soc. Am. B/Vol.9,No.3/March 1992' WRITE(OUTCHAN, *) 'Refractive indices of orthorhombic KNbO3' WRITE(OUTCHAN, *) 'B. Zysset' ELSE IF ( CRYSTAL .EQ. 233) THEN WRITE(OUTCHAN, *) 'Crystal: KNbO3 (100 C )' WRITE(OUTCHAN, *) 'J. Opt. Soc. Am. B/Vol.9,No.3/March 1992' WRITE(OUTCHAN, *) 'Refractive indices of orthorhombic KNbO3' WRITE(OUTCHAN, *) 'B. Zysset' ELSE IF ( CRYSTAL .EQ. 234) THEN WRITE(OUTCHAN, *) 'Crystal: KNbO3 (140 C )' WRITE(OUTCHAN, *) 'J. Opt. Soc. Am. B/Vol.9,No.3/March 1992' WRITE(OUTCHAN, *) 'Refractive indices of orthorhombic KNbO3' WRITE(OUTCHAN, *) 'B. Zysset' ELSE IF ( CRYSTAL .EQ. 235) THEN WRITE(OUTCHAN, *) 'Crystal: KNbO3 (180 C )' WRITE(OUTCHAN, *) 'J. Opt. Soc. Am. B/Vol.9,No.3/March 1992' WRITE(OUTCHAN, *) 'Refractive indices of orthorhombic KNbO3' WRITE(OUTCHAN, *) 'B. Zysset' ELSE WRITE(OUTCHAN, *) 'Crystal: ????' ENDIF LAMBDA=1 CALL INDEX(CRYSTAL, LAMBDA, NX, NY, NZ,ABS_MIN, ABS_MAX) WRITE(OUTCHAN, "(' Lower valid wavelength (microns) = ',F10.4)") ABS_MIN WRITE(OUTCHAN, "(' Upper valid wavelength (microns) = ',F10.4)") ABS_MAX WRITE(*, "(' Lower valid wavelength (microns) = ',F10.4)") ABS_MIN WRITE(*, "(' Upper valid wavelength (microns) = ',F10.4)") ABS_MAX END !*************************************************************************** SUBROUTINE INITFILE(FILENAME, OUTCHAN,INCHAN, CRYSTAL,MENU) IMPLICIT NONE CHARACTER*30 PROGRAMNAME CHARACTER*30 FILENAME ! Output channel INTEGER OUTCHAN,INCHAN INTEGER CRYSTAL,MENU,MENU2 COMMON /CM2/ PROGRAMNAME MENU2 = 0 WRITE(*, *) 'How will the parameters for the calculation be chosen?' WRITE(*, *) "" WRITE(*, *) '(1) Enter parameters manually' WRITE(*, *) "" WRITE(*, *) '(2) Read parameters from text file' WRITE(*, *) ' (see help menu for file format)' WRITE(*, *) "" WRITE(*, *) 'Choice:' READ(*, *) MENU2 IF (MENU2==1) THEN WRITE(*, "(' Output File name =',$)") READ(*, *) FILENAME OUTCHAN = 23 INCHAN = 5 OPEN(UNIT=OUTCHAN, FILE=FILENAME, STATUS='NEW') WRITE(OUTCHAN, "(' Program: ', A30)") PROGRAMNAME WRITE(OUTCHAN, "(' Output file: ', A30)") FILENAME IF (MENU .NE. 4 )THEN CALL LIST WRITE(*, "( ' Crystal number = ' , $)") READ(*, *) CRYSTAL CALL FLIST(OUTCHAN,CRYSTAL) ENDIF ELSE INCHAN = 33 OUTCHAN = 23 IF (MENU .EQ. 1) THEN OPEN(UNIT=INCHAN, FILE='PHASEMATCHIN1', STATUS='OLD') ELSEIF ( MENU .EQ. 2) THEN OPEN(UNIT=INCHAN, FILE='PHASEMATCHIN2', STATUS='OLD') ELSEIF ( MENU .EQ. 3) THEN OPEN(UNIT=INCHAN, FILE='PHASEMATCHIN3', STATUS='OLD') ELSEIF ( MENU .EQ. 4) THEN OPEN(UNIT=INCHAN, FILE='PHASEMATCHIN4', STATUS='OLD') ELSEIF ( MENU .EQ. 5) THEN OPEN(UNIT=INCHAN, FILE='PHASEMATCHIN5', STATUS='OLD') ELSEIF ( MENU .EQ. 6) THEN OPEN(UNIT=INCHAN, FILE='PHASEMATCHIN6', STATUS='OLD') ELSEIF ( MENU .EQ. 7) THEN OPEN(UNIT=INCHAN, FILE='PHASEMATCHIN7', STATUS='OLD') ELSEIF ( MENU .EQ. 8) THEN OPEN(UNIT=INCHAN, FILE='PHASEMATCHIN8', STATUS='OLD') ELSEIF ( MENU .EQ. 9) THEN OPEN(UNIT=INCHAN, FILE='PHASEMATCHIN9', STATUS='OLD') ELSEIF ( MENU .EQ. 10) THEN OPEN(UNIT=INCHAN, FILE='PHASEMATCHIN10', STATUS='OLD') ENDIF READ(INCHAN, *) FILENAME OPEN(UNIT=OUTCHAN, FILE=FILENAME, STATUS='NEW') WRITE(OUTCHAN, "(' Program: ', A30)") PROGRAMNAME WRITE(OUTCHAN, "(' Output file: ', A30)") FILENAME IF (MENU .NE. 4 )THEN CALL LIST READ(INCHAN, *) CRYSTAL CALL FLIST(OUTCHAN,CRYSTAL) ENDIF ENDIF END !*************************************************************************** SUBROUTINE INDEX(CRYSTAL, LAMBDA, NX, NY, NZ, & ABS_MIN, ABS_MAX) IMPLICIT NONE INTEGER CRYSTAL REAL*8 LAMBDA REAL*8 ABS_MIN,ABS_MAX REAL*8 NX,NY,NZ ! *************************UNIAXIAL******************************** ! Ag3AsS3 ( proustite ) ! Rivista Del Nuovo Cimento Vol 6, N 3, p 295-319 (1976 ) ! " Nonlinear Materials for Tunable Coherent Light Sources ! in the Middle Infra-Red " ! G. AIROLDI IF ( CRYSTAL .EQ. 10 ) THEN NX = SQRT( 9.220 + .4454 / ( LAMBDA**2 - 0.1264 ) & + 1733. / (LAMBDA**2 - 1000.) ) NY = NX NZ = SQRT( 7.007 + .3230 / ( LAMBDA**2 - 0.1192 ) & + 660. / ( LAMBDA**2 - 1000.) ) ABS_MIN = 0.5 ABS_MAX = 10. ! AgGaSe2 Ref 1 ! Kildal, H., J. Mikkelsen, Opt. Commun., 9, 315, (1973) ! ~0.8 - 17 (35% Refl. loss) ! 1 - 13.5 ELSE IF ( CRYSTAL .EQ. 20 ) THEN NX = SQRT( 3.9362 + 2.9113/(1 - (0.38821/LAMBDA)**2) + & 1.7954/(1 - (40/LAMBDA)**2) ) NZ = SQRT( 3.3132 + 3.3616/(1 - (0.38201/LAMBDA)**2) + & 1.7677/(1 - (40/LAMBDA)**2) ) NY = NX ABS_MIN = 0.8 ABS_MAX = 17.0 ! AgGaSe2 Ref 2 ! Bhar, G. C., Appl. Opt., 15, 305 (1976) ! ~0.8 - 17 (35% Refl. loss) ! 0.7 - 13.5 ELSE IF ( CRYSTAL .EQ. 30 ) THEN NX = SQRT( 4.6453 + 2.2057/(1 - (0.43347/LAMBDA)**2) + & 1.8377/(1-(40/LAMBDA)**2) ) NZ = SQRT( 5.2912 + 1.3970/(1 - (0.53339/LAMBDA)**2) + & 1.9282/(1 - (40/LAMBDA)**2) ) NY = NX ABS_MIN = 0.8 ABS_MAX = 17.0 ! BBO Ref 1 ! Katz, K., IEEE J. Quant. Electron., 22, 1013-4 (1986) ! ~0.2 - 2.1 ! 0.3 - 1.5 ? ELSE IF ( CRYSTAL .EQ. 40 ) THEN NX = SQRT( 2.7359 + 0.01878/(LAMBDA**2 - 0.01822) - & 0.01354*LAMBDA**2 ) NZ = SQRT( 2.3753 + 0.01224/(LAMBDA**2 - 0.01667) - & 0.01516*LAMBDA**2 ) NY = NX ABS_MIN = 0.2 ABS_MAX = 2.1 ! BBO Ref 2 ! Graham, A. Zalkin, J. Appl. Phys., 62 (1987) ! 0.2-2.8 ! 0.4 - 1.0 ELSE IF ( CRYSTAL .EQ. 50 ) THEN IF ( LAMBDA .LE. 12.0 ) THEN NX = SQRT( 2.7405 + 0.0184/(LAMBDA**2 - 0.0179) - & 0.0155*LAMBDA**2 ) NZ = SQRT( 2.3730 + 0.0128/(LAMBDA**2 - 0.0156) - & 0.0044*LAMBDA**2 ) ELSE NX = 0.2 NZ = 2.8 ENDIF NY = NX ABS_MIN = 0.3 ABS_MAX = 5.0 ! GaSe ! Abdullaev, G.B., L.A. Kulevskii, A.M. Prokharov, A.D. Savel'ev, ! E. Yu. Salaev, V.V. Smirnov, JETP Lett., 16, 90 (1972) ! ~0.65 - 18 (40% Relf. loss) ! 0.63 - 10.6 ELSE IF ( CRYSTAL .EQ. 60 ) THEN NX = SQRT( -0.05466/LAMBDA**4 + 0.48605/LAMBDA**2 + & 7.8902 - (0.000824)*LAMBDA**2 - & (0.00000273)*LAMBDA**4 ) NZ = SQRT( 6.0476 + 0.3423/(LAMBDA**2 - 0.16491) - & (0.001042)*LAMBDA**2 ) NY = NX ABS_MIN = 0.65 ABS_MAX = 18.0 ! KDP Ref 1 ! Zernike, Frits, Jr., J. Opt. Soc. Am., 54, 1215 (1964) ! ~0.2 - 1.5 ! 0.21 - 1.5 ELSE IF ( CRYSTAL .EQ. 70 ) THEN NX = SQRT( 0.01011279/(LAMBDA**2 - 0.01294238) + & 0.03249268/(0.0025 - 1/LAMBDA**2) + 2.260476 ) NZ = SQRT( 0.00865325/(LAMBDA**2 - 0.01229326) + & 0.00806984/(0.0025 - 1/LAMBDA**2) + 2.133831 ) NY = NX ABS_MIN = 0.2 ABS_MAX = 1.5 ! KDP Ref 2 ! Barnes, N.P., D.J. Gettemy, R.S. Adhav, J. Opt. Soc. Am., ! 72, 895 (1982) ! ~0.2 - 1.5 ! 0.41 - 0.65 ELSE IF ( CRYSTAL .EQ. 80 ) THEN NX = SQRT( 1.0 + 1.24361*LAMBDA**2/(LAMBDA**2 - 0.00959) ) NZ = SQRT( 1.0 + 1.12854*LAMBDA**2/(LAMBDA**2 - 0.00841) ) NY = NX ABS_MIN = 0.2 ABS_MAX = 1.5 ! KD*P Ref 1 ! Barnes, N.P., D.J. Gettemy, R.S. Adhav, J. Opt. Soc. Am., ! 72, 895 (1982) ! ~0.2 - 1.8 (98% composition) ! 0.41 - 0.65 ELSE IF ( CRYSTAL .EQ. 90 ) THEN NX = SQRT( 1.0 + 1.23318*LAMBDA**2/(LAMBDA**2 - 0.00931) ) NZ = SQRT( 1.0 + 1.12354*LAMBDA**2/(LAMBDA**2 - 0.00828) ) NY = NX ABS_MIN = 0.2 ABS_MAX = 1.8 ! KD*P Ref 2 ! Phillips, R.A., J. Opt. Soc. Am., 56, 629 (1966) ! ~0.2 - 1.8 (98% composition) ! ? ELSE IF ( CRYSTAL .EQ. 100 ) THEN NX = SQRT( 1.661824 + 0.585337/(1 - (0.016017/LAMBDA**2)) + & 0.691221/(1 - (30.0/LAMBDA**2)) ) NZ = SQRT( 1.687522 + 0.447488/(1 - (0.017039/LAMBDA**2)) + & 0.596216/(1 - (30.0/LAMBDA**2)) ) NY = NX ABS_MIN = 0.2 ABS_MAX = 1.8 ! K3Li2Nb5O15 ! Nikogosyan, Handbook of Nonlinear Optical Crystals, 1991, p. 76 ! 0.4 - 5 ! 0.45 - 1? ELSE IF ( CRYSTAL .EQ. 110 ) THEN NX = SQRT( 1.0 + 3.708*LAMBDA**2/(LAMBDA**2 - 0.04601) ) NZ = SQRT( 1.0 + 3.349*LAMBDA**2/(LAMBDA**2 - 0.03564) ) NY = NX ABS_MIN = 0.4 ABS_MAX = 5.0 ! LiIO3 Ref 6 (use Ref. 7) ! Takizawa, K., M. Okada, S. Ieiri, Opt. Commun., 23, 279 (1977) ! 0.3 - > 5.0 ! 0.35 - 0.65 ELSE IF ( CRYSTAL .EQ. 120 ) THEN NX = SQRT( 3.4095 + 0.047664/(LAMBDA**2 - 0.033991) ) NZ = SQRT( 2.9163 + 0.034514/(LAMBDA**2 - 0.031034) ) NY = NX ABS_MIN = 0.3 ABS_MAX = 5.0 ! LiIO3 Ref 7 ! Phillips, R.A., J. Opt. Soc. Am., 56, 629 (1966) ! 0.3 - > 5.0 ! 0.35 - 5.0 ELSE IF ( CRYSTAL .EQ. 130 ) THEN NX = SQRT( 2.03132 + 1.37623/(1 - (0.0350832/LAMBDA**2)) + & 1.06745/(1 - (169/LAMBDA**2)) ) NZ = SQRT( 1.83086 + 1.08807/(1 - (0.031381/LAMBDA**2)) + & 0.554582/(1 - (158.76/LAMBDA**2)) ) NY = NX ABS_MIN = 0.3 ABS_MAX = 5.0 ! LiNO3 ! Smith, D.S., H.D. Riccius, R.P. Edwin, Opt. Commun., 17, ! 332, (1976) ! Probably to 3.4 ! 0.4 - 3.39 ELSE IF ( CRYSTAL .EQ. 140 ) THEN NX = SQRT( 4.9048 - 0.11768/(0.04750 - LAMBDA**2) - & 0.027169*LAMBDA**2 ) NZ = SQRT( 4.5820 - 0.099169/(0.044432 - LAMBDA**2) - & 0.021950*LAMBDA**2 ) NY = NX ABS_MIN = 0.4 ABS_MAX = 3.4 ! Tl3AsSe3 ! Feichtner, J.D., G.W. Roland, Appl.Opt., 11, 993 (1972) ! 1.26 - 16 (45% not corrected for refl.) ! 1.5 - 5.3 ELSE IF ( CRYSTAL .EQ. 150 ) THEN NX = SQRT( 1.0 + 9.977*LAMBDA**2/(LAMBDA**2 - (0.435)**2) + & 0.067*LAMBDA**2/(LAMBDA**2 - (20)**2) ) NZ = SQRT( 1.0 + 8.783*LAMBDA**2/(LAMBDA**2 - (0.435)**2) + & 0.051*LAMBDA**2/(LAMBDA**2 - (20)**2) ) NY = NX ABS_MIN = 1.26 ABS_MAX = 16.0 ! ZnGeP2 ! Dmitriev p. 85 ! 0.74 - 12 ELSE IF ( CRYSTAL .EQ. 160 ) THEN NX = SQRT( 4.4733 + 5.26575*LAMBDA**2/(LAMBDA**2 - 0.13381) + & 1.49085*LAMBDA**2/(LAMBDA**2 - 662.55 ) ) NZ = SQRT( 4.63318 + 5.34215*LAMBDA**2/(LAMBDA**2 - 0.142555) + & 1.45785*LAMBDA**2/(LAMBDA**2 - 662.55) ) NY = NX ABS_MIN = 0.74 ABS_MAX = 12.0 ! AgGaS2 ! G. C. Bahr, Appl. Opt. 15, 305 (1976) ! 0.5 - 13 ELSE IF ( CRYSTAL .EQ. 170 ) THEN NX = SQRT( 3.6280 + 2.1686*LAMBDA**2/(LAMBDA**2 - 0.1003) + & 2.1753*LAMBDA**2/(LAMBDA**2 - 950 ) ) NZ = SQRT( 4.0172+ 1.5274*LAMBDA**2/(LAMBDA**2 - 0.1310) + & 2.1699*LAMBDA**2/(LAMBDA**2 - 950) ) NY = NX ABS_MIN = 0.5 ABS_MAX = 13.0 !***************************BIAXIAL************************************** ! KTP ! H. Vanherzeele, J. D. Bierlein, F. C. Zumsteg, Appl. Opt., 27, ! 3314 (1988) ELSE IF ( CRYSTAL .EQ. 200 ) THEN NX = SQRT( 2.1146 + 0.89188/(1 - (0.20861/LAMBDA)**2) - (0.01320*LAMBDA**2)) NY = SQRT( 2.1518 + 0.87862/(1 - (0.21801/LAMBDA)**2) - (0.01327*LAMBDA**2)) NZ = SQRT( 2.3136 + 1.00012/(1 - (0.23831/LAMBDA)**2) - (0.01679*LAMBDA**2)) ABS_MIN = 0.35 ABS_MAX = 4.50 ! KTA ! D. L. Fenimore, K. L. Schepler, U. B. Ramabadran, ! S. R. McPherson, "Infrared corrected Sellmeier coefficients ! for potassium titanyl arsenate," J. Opt. Soc. B, 12 (5), ! 794-796 (May 1995) ELSE IF ( CRYSTAL .EQ. 210 ) THEN NX = SQRT( 1.90713 + 1.23522/(1 - (0.19692/LAMBDA)**2) - (0.01025*LAMBDA**2)) NY = SQRT( 2.15912 + 1.00099/(1 - (0.21844/LAMBDA)**2) - (0.01096*LAMBDA**2)) NZ = SQRT( 2.14786 + 1.29559/(1 - (0.22719/LAMBDA)**2) - (0.01436*LAMBDA**2)) ABS_MIN = 0.35 ABS_MAX = 4.0 ! Magnesium-Barium Flouride ! V. G. Dmitriev, G. G. Gurzadyan, D. N. Nikogosyan, ! "Handbook of Nonlinear Crystals," (Springer-Verlag, ! Berlin) 1991, pp.112-113 ELSE IF ( CRYSTAL .EQ. 220 ) THEN NX = SQRT( 2.077 + 0.0076/(LAMBDA**2 - 0.0079) ) NY = SQRT( 2.1238 + 0.0086/(LAMBDA**2) ) NZ = SQRT( 2.1462 + 0.00736/(LAMBDA**2 - 0.009) ) ABS_MIN = 0.3 ABS_MAX = 5.0 ! KNbO3 ( 22ūC ) ! J. Opt. Soc. Am. B/Vol.9,No.3/March 1992 ! Refractive indices of orthorhombic KNbO3 ! B. Zysset ELSE IF ( CRYSTAL .EQ. 230 ) THEN NX = SQRT( 1. + & 20.05519 * LAMBDA**2 * & 6.6645268D-2 /( LAMBDA**2 - & 6.6645268D-2 )+ 149.8408* & LAMBDA**2 *1.6664790D-2 / & ( LAMBDA**2 - & 1.6664790D-2 )- & 2.517432D-2 * LAMBDA**2 ) NY = SQRT( 1. + & 19.37347 * LAMBDA**2 * & 7.4390211D-2 /( LAMBDA**2 - & 7.4390211D-2 )+ & 135.4992 * LAMBDA**2 * & 1.8770074D-2 /( LAMBDA**2 - & 1.8770074D-2 )- & 2.845018D-2 * LAMBDA**2 ) NZ = SQRT( 1. + & 16.0917 * LAMBDA**2 * & 6.5141246D-2 /( LAMBDA**2 - & 6.5141246D-2 )+ & 165.4431 * LAMBDA**2 * & 1.4331497D-2 /( LAMBDA**2 - & 1.4331497D-2 )- & 1.943289D-2 * LAMBDA**2 ) ABS_MIN = .4 ABS_MAX = 4.5 ! KNbO3 ( 50ūC ) ! J. Opt. Soc. Am. B/Vol.9,No.3/March 1992 ! Refractive indices of orthorhombic KNbO3 ! B. Zysset ELSE IF ( CRYSTAL .EQ. 231 ) THEN NX = SQRT( 1. + & 20.27475 * LAMBDA**2 * & 6.6907198D-2 /( LAMBDA**2 - & 6.6907198D-2 )+ & 152.338 * LAMBDA**2 * & 1.6279778D-2 /( LAMBDA**2 - & 1.6279778D-2 )- & 2.514857D-2 * LAMBDA**2 ) NY = SQRT( 1. + & 20.6625 * LAMBDA**2 * & 7.3657362D-2 /( LAMBDA**2 - & 7.3657362D-2 )+ & 144.1607 * LAMBDA**2 * & 1.7054649D-2 /( LAMBDA**2 - & 1.7054649D-2 )- & 2.827127D-2 * LAMBDA**2 ) NZ = SQRT( 1. + & 18.24447 * LAMBDA**2 * & 2.5275264D-1 /( LAMBDA**2 - & 2.5275264D-1 )+ & 181.1303 * LAMBDA**2 * & 1.2486341D-2 /( LAMBDA**2 - & 1.2486341D-2 )- & 1.945428D-2 * LAMBDA**2 ) ABS_MIN = .4 ABS_MAX = 4.5 ! KNbO3 ( 75ūC ) ! J. Opt. Soc. Am. B/Vol.9,No.3/March 1992 ! Refractive indices of orthorhombic KNbO3 ! B. Zysset ELSE IF ( CRYSTAL .EQ. 232 ) THEN NX = SQRT( 1. + & 19.33691 * LAMBDA**2 * & 6.8002224D-2 /( LAMBDA**2 - & 6.8002224D-2 )+ & 150.4971 * LAMBDA**2 * & 1.6771525D-2 /( LAMBDA**2 - & 1.6771525D-2 )- & 2.514893D-2 * LAMBDA**2 ) NY = SQRT( 1. + & 20.06891 * LAMBDA**2 * & 7.4272983D-2 /( LAMBDA**2 - & 7.4272983D-2 )+ & 142.7868 * LAMBDA**2 * & 1.7410769D-2 /( LAMBDA**2 - & 1.7410769D-2 )- & 2.80751D-2 * LAMBDA**2 ) NZ = SQRT( 1. + & 17.56954 * LAMBDA**2 * & 6.4859835D-2 /( LAMBDA**2 - & 6.4859835D-2 )+ & 176.6712 * LAMBDA**2 * & 1.2989853D-2 /( LAMBDA**2 - & 1.2989853D-2 )- & & 1.95288D-2 * LAMBDA**2 ) ABS_MIN = .4 ABS_MAX = 4.5 ! KNbO3 (100ūC ) ! J. Opt. Soc. Am. B/Vol.9,No.3/March 1992 ! Refractive indices of orthorhombic KNbO3 ! B. Zysset ELSE IF ( CRYSTAL .EQ. 233 ) THEN NX = SQRT( 1. + & 16.60832 * LAMBDA**2 * & 7.0558671D-2 /( LAMBDA**2 - & 7.0558671D-2 )+ & 142.0508 * LAMBDA**2 * & 1.8797270D-2 /( LAMBDA**2 - & 1.8797270D-2 )- & 2.51586D-2 * LAMBDA**2 ) NY = SQRT( 1. + & 14.37713 * LAMBDA**2 * & 7.9499205D-2 /( LAMBDA**2 - & 7.9499205D-2 )+ & 123.3481 * LAMBDA**2 * & 2.2931418D-2 /( LAMBDA**2 - & 2.2931418D-2 )- & 2.79633D-2 * LAMBDA**2 ) NZ = SQRT( 1. + & 14.0194 * LAMBDA**2 * & 6.8759235D-2 /( LAMBDA**2 - & 6.8759235D-2 )+ & 158.4156 * LAMBDA**2 * & 1.5650029D-2 /( LAMBDA**2 - & 1.5650029D-2 )- & 1.968455D-2 * LAMBDA**2 ) ABS_MIN = .4 ABS_MAX = 4.5 ! KNbO3 (140ūC ) ! J. Opt. Soc. Am. B/Vol.9,No.3/March 1992 ! Refractive indices of orthorhombic KNbO3 ! B. Zysset ELSE IF ( CRYSTAL .EQ. 234 ) THEN NX = SQRT( 1. + 16.04692 * LAMBDA**2 * & 7.1705899D-2 /( LAMBDA**2 - 7.1705899D-2 )+ & 140.7942 * LAMBDA**2 * 1.9153316D-2 /( LAMBDA**2 - & 1.9153316D-2 )- 2.508202D-2 * LAMBDA**2 ) NY = SQRT( 1. + & 18.87445 * LAMBDA**2 * & 7.5550311D-2 /( LAMBDA**2 - & 7.5550311D-2 )+ & 138.8675 * LAMBDA**2 * & 1.8301250D-2 /( LAMBDA**2 - & 1.8301250D-2 )- & 2.771253D-2 * LAMBDA**2 ) NZ = SQRT( 1. + & 16.48284 * LAMBDA**2 * & 6.6932302D-2 /( LAMBDA**2 - & 6.6932302D-2 )+ & 167.7579 * LAMBDA**2 * & 1.4040169D-2 /( LAMBDA**2 - & 1.4040169D-2 )- & 1.982648D-2 * LAMBDA**2 ) ABS_MIN = .4 ABS_MAX = 4.5 ! KNbO3 (180ūC ) ! J. Opt. Soc. Am. B/Vol.9,No.3/March 1992 ! Refractive indices of orthorhombic KNbO3 ! B. Zysset ELSE IF ( CRYSTAL .EQ. 235 ) THEN NX = SQRT( 1. + & 13.89662 * LAMBDA**2 * & 7.4910738D-2 /( LAMBDA**2 - & 7.4910738D-2 )+ & 137.2157 * LAMBDA**2 * & 2.0502142D-2 /( LAMBDA**2 - & 2.0502142D-2 )- & 2.511282D-2 * LAMBDA**2 ) NY = SQRT( 1. + & 22.1202 * LAMBDA**2 * & 7.3449810D-2 /( LAMBDA**2 - & 7.3449810D-2 )+ & 162.8852 * LAMBDA**2 * & 1.4350390D-2 /( LAMBDA**2 - & 1.4350390D-2 )- & 2.762299D-2 * LAMBDA**2 ) NZ = SQRT( 1. + & 21.47866 * LAMBDA**2 * & 6.3894396D-2 /( LAMBDA**2 - & 6.3894396D-2 )+ & 215.6493 * LAMBDA**2 * & 9.7616759D-3 /( LAMBDA**2 - & 9.7616759D-3 )- & 2.004523D-2 * LAMBDA**2 ) ABS_MIN = .4 ABS_MAX = 4.5 ELSE CALL ERR_MSG(" Requested crystal is not available.", 45) ENDIF END SUBROUTINE UNCMND (N,X0,FCN,X,F,INFO,W,LW) !***BEGIN PROLOGUE UNCMND !***DATE WRITTEN 870923 (YYMMDD) !***REVISION DATE 871222 (YYMMDD) !***CATEGORY NO. G1B1A1 !***KEYWORDS UNCONSTRAINED MINIMIZATION !***AUTHOR NASH, S.G., (GEORGE MASON UNIVERSITY) !***PURPOSE UNCMND minimizes a smooth nonlinear function of n variables. ! A subroutine that computes the function value at any point ! must be supplied, but derivative values are not required. ! UNCMND provides a simple interface to more flexible lower ! level routines. User has no control over options. ! !***DESCRIPTION ! From the book, "Numerical Methods and Software" by ! D. Kahaner, C. Moler, S. Nash ! Prentice Hall, 1988 ! ! This routine uses a quasi-Newton algorithm with line search ! to minimize the function represented by the subroutine FCN. ! At each iteration, the nonlinear function is approximated ! by a quadratic function derived from a Taylor series. ! The quadratic function is minimized to obtain a search direction, ! and an approximate minimum of the nonlinear function along ! the search direction is found using a line search. The ! algorithm computes an approximation to the second derivative ! matrix of the nonlinear function using quasi-Newton techniques. ! ! The UNCMND package is quite general, and provides many options ! for the user. However, this subroutine is designed to be ! easy to use, with few choices allowed. For example: ! ! 1. Only function values need be computed. First derivative ! values are obtained by finite-differencing. This can be ! very costly when the number of variables is large. ! ! 2. It is assumed that the function values can be obtained ! accurately (to an accuracy comparable to the precision of ! the computer arithmetic). ! ! 3. At most 150 iterations are allowed. ! ! 4. It is assumed that the function values are well-scaled, ! that is, that the optimal function value is not pathologically ! large or small. ! ! For more information, see the reference listed below. ! ! PARAMETERS ! ---------- ! N --> INTEGER ! Dimension of problem ! X0(N) --> DOUBLE PRECISION ! Initial estimate of minimum ! FCN --> Name of routine to evaluate minimization function. ! Must be declared EXTERNAL in calling routine, and ! have calling sequence ! SUBROUTINE FCN(N, X, F) ! with N and X as here, F the computed function value. ! X(N) <-- DOUBLE PRECISION ! Local minimum ! F <-- DOUBLE PRECISION ! Function value at local minimum X ! INFO <-- INTEGER ! Termination code ! INFO = 0: Optimal solution found ! INFO = 1: Terminated with gradient small, ! X is probably optimal ! INFO = 2: Terminated with stepsize small, ! X is probably optimal ! INFO = 3: Lower point cannot be found, ! X is probably optimal ! INFO = 4: Iteration limit (150) exceeded ! INFO = 5: Too many large steps, ! function may be unbounded ! INFO = -1: Insufficient workspace ! W(LW) --> DOUBLE PRECISION ! Workspace ! LW --> INTEGER ! Size of workspace, at least N*(N+10) ! !***REFERENCES R.B. SCHNABEL, J.E. KOONTZ, AND BE.E. WEISS, A MODULAR ! SYSTEM OF ALGORITHMS FOR UNCONSTRAINED MINIMIZATION, ! REPORT CU-CS-240-82, COMP. SCI. DEPT., UNIV. OF ! COLORADO AT BOULDER, 1982. !***ROUTINES CALLED OPTDRD, XERROR !***END PROLOGUE UNCMND IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X0(N),X(N),W(LW) CHARACTER ERRMSG*68 EXTERNAL FCN, D1FCND, D2FCND !---------------------------------------------------------------- ! SUBDIVIDE WORKSPACE !---------------------------------------------------------------- !***FIRST EXECUTABLE STATEMENT UNCMND IG = 1 IT = IG + N IW1 = IT + N IW2 = IW1 + N IW3 = IW2 + N IW4 = IW3 + N IW5 = IW4 + N IW6 = IW5 + N IW7 = IW6 + N IW8 = IW7 + N IA = IW8 + N LWMIN = IA + N*N-1 IF (LWMIN .GT. LW) THEN INFO = -1 WRITE(ERRMSG, '( & ''UNCMND ERROR (INFO=-1) -- INSUFFICIENT WORKSPACE'',& '', LW = '', I5 )' ) LW CALL XERROR(ERRMSG(1:60), 60, -1, 0) RETURN END IF !---------------------------------------------------------------- ! SET UP PARAMETERS FOR OPTDRD !---------------------------------------------------------------- ! PARAMETERS THAT SHOULD NOT BE CHANGED WHEN USING CONDENSED CODE ! ! NR = PARAMETER USED TO DIVIDE WORKSPACE ! METHOD = 1 (LINE SEARCH) -- DO NOT CHANGE ! MSG = 9 => NO PRINTING, N=1 ALLOWED ! IAGFLG = 1 => ANALYTIC GRADIENT SUPPLIED (0 OTHERWISE) ! IAHFLG = 1 => ANALYTIC HESSIAN SUPPLIED (0 OTHERWISE) ! IPR = DEVICE FOR OUTPUT (IRRELEVANT IN CURRENT VERSION) ! DLT = (IRRELEVANT PARAMETER FOR METHOD = 1) ! EPSM = MACHINE EPSILON ! IEXP = 1 => FUNCTION EXPENSIVE TO EVALUATE (IEXP = 0 => CHEAP) ! NR = N METHOD = 1 MSG = 9 IAGFLG = 0 IAHFLG = 0 IPR = 0 DLT = -1.0D0 EPSM = D1MACH(4) IEXP = 1 ! ! PARAMETERS THAT MAY BE CHANGED: ! ! NDIGIT = -1 => OPTDRD ASSUMES F IS FULLY ACCURATE ! ITNLIM = 150 = MAXIMUM NUMBER OF ITERATIONS ALLOWED ! GRADTL = ZERO TOLERANCE FOR GRADIENT, FOR CONVERGENCE TESTS ! STEPMX = MAXIMUM ALLOWABLE STEP SIZE ! STEPTL = ZERO TOLERANCE FOR STEP, FOR CONVERGENCE TESTS ! FSCALE = TYPICAL ORDER-OF-MAGNITUDE SIZE OF FUNCTION ! TYPSIZ = TYPICAL ORDER-OF-MAGNITUDE SIZE OF X (STORED IN W(LT)) ! NDIGIT = -1 ITNLIM = 1500 GRADTL = EPSM**(1.0D0/3.0D0) STEPMX = 0.0D0 STEPTL = SQRT(EPSM) FSCALE = 0.01D0 DO 10 LT = IT,IT+N-1 W(LT) = 0.01D0 10 CONTINUE ! ! MINIMIZE FUNCTION ! CALL OPTDRD (NR, N, X0, FCN, D1FCND, D2FCND, W(IT), FSCALE, & METHOD, IEXP, MSG, NDIGIT, ITNLIM, IAGFLG, IAHFLG, & IPR, DLT, GRADTL, STEPMX, STEPTL, & X, F, W(IG), INFO, W(IA), & W(IW1), W(IW2), W(IW3), W(IW4), & W(IW5), W(IW6), W(IW7), W(IW8)) ! IF (INFO .EQ. 1) THEN WRITE(ERRMSG, '( & ''UNCMND WARNING -- INFO = 1'', & '': PROBABLY CONVERGED, GRADIENT SMALL'')' ) CALL XERROR(ERRMSG(1:62), 62, INFO, 0) END IF IF (INFO .EQ. 2) THEN WRITE(ERRMSG, '( & ''UNCMND WARNING -- INFO = 2'', & '': PROBABLY CONVERGED, STEPSIZE SMALL'')' ) CALL XERROR(ERRMSG(1:62), 62, INFO, 0) END IF IF (INFO .EQ. 3) THEN WRITE(ERRMSG, '( & ''UNCMND WARNING -- INFO = 3'', & '': CANNOT FIND LOWER POINT'')' ) CALL XERROR(ERRMSG(1:51), 51, INFO, 0) END IF IF (INFO .EQ. 4) THEN WRITE(ERRMSG, '( & ''UNCMND WARNING -- INFO = 4'', & '': TOO MANY ITERATIONS'')' ) CALL XERROR(ERRMSG(1:47), 47, INFO, 0) END IF IF (INFO .EQ. 5) THEN WRITE(ERRMSG, '( & ''UNCMND WARNING -- INFO = 5'', & '': TOO MANY LARGE STEPS, POSSIBLY UNBOUNDED'')' ) CALL XERROR(ERRMSG(1:68), 68, INFO, 0) END IF ! RETURN END SUBROUTINE XERROR(MESSG,NMESSG,NERR,LEVEL) !***BEGIN PROLOGUE XERROR !***DATE WRITTEN 790801 (YYMMDD) !***REVISION DATE 820801 (YYMMDD) !***CATEGORY NO. R3C !***KEYWORDS ERROR,XERROR PACKAGE !***AUTHOR JONES, R. E., (SNLA) !***PURPOSE Processes an error (diagnostic) message. !***DESCRIPTION ! From the book "Numerical Methods and Software" ! by D. Kahaner, C. Moler, S. Nash ! Prentice Hall 1988 ! Abstract ! XERROR processes a diagnostic message, in a manner ! determined by the value of LEVEL and the current value ! of the library error control flag, KONTRL. ! (See subroutine XSETF for details.) ! ! Description of Parameters ! --Input-- ! MESSG - the Hollerith message to be processed, containing ! no more than 72 characters. ! NMESSG- the actual number of characters in MESSG. ! NERR - the error number associated with this message. ! NERR must not be zero. ! LEVEL - error category. ! =2 means this is an unconditionally fatal error. ! =1 means this is a recoverable error. (I.e., it is ! non-fatal if XSETF has been appropriately called.) ! =0 means this is a warning message only. ! =-1 means this is a warning message which is to be ! printed at most once, regardless of how many ! times this call is executed. ! ! Examples ! CALL XERROR('SMOOTH -- NUM WAS ZERO.',23,1,2) ! CALL XERROR('INTEG -- LESS THAN FULL ACCURACY ACHIEVED.', ! 43,2,1) ! CALL XERROR('ROOTER -- ACTUAL ZERO OF F FOUND BEFORE INTERVAL F ! 1ULLY COLLAPSED.',65,3,0) ! CALL XERROR('EXP -- UNDERFLOWS BEING SET TO ZERO.',39,1,-1) ! ! Latest revision --- 19 MAR 1980 ! Written by Ron Jones, with SLATEC Common Math Library Subcommittee !***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- ! HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, ! 1982. !***ROUTINES CALLED XERRWV !***END PROLOGUE XERROR CHARACTER*(*) MESSG !***FIRST EXECUTABLE STATEMENT XERROR CALL XERRWV(MESSG,NMESSG,NERR,LEVEL,0,0,0,0,0.,0.) RETURN END SUBROUTINE XERRWV(MESSG,NMESSG,NERR,LEVEL,NI,I1,I2,NR,R1,R2) !***BEGIN PROLOGUE XERRWV !***DATE WRITTEN 800319 (YYMMDD) !***REVISION DATE 870916 (YYMMDD) !***CATEGORY NO. R3C !***KEYWORDS ERROR,XERROR PACKAGE !***AUTHOR JONES, R. E., (SNLA) !***PURPOSE Processes error message allowing 2 integer and two real ! values to be included in the message. !***DESCRIPTION ! From the book "Numerical Methods and Software" ! by D. Kahaner, C. Moler, S. Nash ! Prentice Hall 1988 ! Abstract ! XERRWV processes a diagnostic message, in a manner ! determined by the value of LEVEL and the current value ! of the library error control flag, KONTRL. ! (See subroutine XSETF for details.) ! In addition, up to two integer values and two real ! values may be printed along with the message. ! ! Description of Parameters ! --Input-- ! MESSG - the Hollerith message to be processed. ! NMESSG- the actual number of characters in MESSG. ! NERR - the error number associated with this message. ! NERR must not be zero. ! LEVEL - error category. ! =2 means this is an unconditionally fatal error. ! =1 means this is a recoverable error. (I.e., it is ! non-fatal if XSETF has been appropriately called.) ! =0 means this is a warning message only. ! =-1 means this is a warning message which is to be ! printed at most once, regardless of how many ! times this call is executed. ! NI - number of integer values to be printed. (0 to 2) ! I1 - first integer value. ! I2 - second integer value. ! NR - number of real values to be printed. (0 to 2) ! R1 - first real value. ! R2 - second real value. ! ! Examples ! CALL XERRWV('SMOOTH -- NUM (=I1) WAS ZERO.',29,1,2, ! 1 1,NUM,0,0,0.,0.) ! CALL XERRWV('QUADXY -- REQUESTED ERROR (R1) LESS THAN MINIMUM ( ! 1R2).,54,77,1,0,0,0,2,ERRREQ,ERRMIN) ! ! Latest revision --- 16 SEPT 1987 ! Written by Ron Jones, with SLATEC Common Math Library Subcommittee !***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- ! HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, ! 1982. !***ROUTINES CALLED FDUMP,I1MACH,J4SAVE,XERABT,XERCTL,XERPRT,XERSAV, ! XGETUA !***END PROLOGUE XERRWV CHARACTER*(*) MESSG CHARACTER*20 LFIRST CHARACTER*37 FORM DIMENSION LUN(5) ! GET FLAGS !***FIRST EXECUTABLE STATEMENT XERRWV LKNTRL = J4SAVE(2,0,.FALSE.) MAXMES = J4SAVE(4,0,.FALSE.) ! CHECK FOR VALID INPUT IF ((NMESSG.GT.0).AND.(NERR.NE.0).AND. & (LEVEL.GE.(-1)).AND.(LEVEL.LE.2)) GO TO 10 IF (LKNTRL.GT.0) CALL XERPRT('FATAL ERROR IN...',17) CALL XERPRT('XERROR -- INVALID INPUT',23) IF (LKNTRL.GT.0) CALL FDUMP IF (LKNTRL.GT.0) CALL XERPRT('JOB ABORT DUE TO FATAL ERROR.', & 29) IF (LKNTRL.GT.0) CALL XERSAV(' ',0,0,0,KDUMMY) CALL XERABT('XERROR -- INVALID INPUT',23) RETURN 10 CONTINUE ! RECORD MESSAGE JUNK = J4SAVE(1,NERR,.TRUE.) CALL XERSAV(MESSG,NMESSG,NERR,LEVEL,KOUNT) ! LET USER OVERRIDE LFIRST = MESSG LMESSG = NMESSG LERR = NERR LLEVEL = LEVEL CALL XERCTL(LFIRST,LMESSG,LERR,LLEVEL,LKNTRL) ! RESET TO ORIGINAL VALUES LMESSG = NMESSG LERR = NERR LLEVEL = LEVEL LKNTRL = MAX0(-2,MIN0(2,LKNTRL)) MKNTRL = IABS(LKNTRL) ! DECIDE WHETHER TO PRINT MESSAGE IF ((LLEVEL.LT.2).AND.(LKNTRL.EQ.0)) GO TO 100 IF (((LLEVEL.EQ.(-1)).AND.(KOUNT.GT.MIN0(1,MAXMES))) & .OR.((LLEVEL.EQ.0) .AND.(KOUNT.GT.MAXMES)) & .OR.((LLEVEL.EQ.1) .AND.(KOUNT.GT.MAXMES).AND.(MKNTRL.EQ.1)) & .OR.((LLEVEL.EQ.2) .AND.(KOUNT.GT.MAX0(1,MAXMES)))) GO TO 100 IF (LKNTRL.LE.0) GO TO 20 CALL XERPRT(' ',1) ! INTRODUCTION IF (LLEVEL.EQ.(-1)) CALL XERPRT & ('WARNING MESSAGE...THIS MESSAGE WILL ONLY BE PRINTED ONCE.',57) IF (LLEVEL.EQ.0) CALL XERPRT('WARNING IN...',13) IF (LLEVEL.EQ.1) CALL XERPRT & ('RECOVERABLE ERROR IN...',23) IF (LLEVEL.EQ.2) CALL XERPRT('FATAL ERROR IN...',17) 20 CONTINUE ! MESSAGE CALL XERPRT(MESSG,LMESSG) CALL XGETUA(LUN,NUNIT) ISIZEI = LOG10(FLOAT(I1MACH(9))) + 1.0 ISIZEF = LOG10(FLOAT(I1MACH(10))**I1MACH(11)) + 1.0 DO 50 KUNIT=1,NUNIT IUNIT = LUN(KUNIT) ! IF (IUNIT.EQ.0) IUNIT = I1MACH(4) DO 22 I=1,MIN(NI,2) WRITE (FORM,21) I,ISIZEI 21 FORMAT ('(11X,21HIN ABOVE MESSAGE, I',I1,'=,I',I2,') ') IF(IUNIT.EQ.0)THEN IF (I.EQ.1) WRITE (*,FORM) I1 IF (I.EQ.2) WRITE (*,FORM) I2 ELSE IF (I.EQ.1) WRITE (IUNIT,FORM) I1 IF (I.EQ.2) WRITE (IUNIT,FORM) I2 ENDIF 22 CONTINUE DO 24 I=1,MIN(NR,2) WRITE (FORM,23) I,ISIZEF+10,ISIZEF 23 FORMAT ('(11X,21HIN ABOVE MESSAGE, R',I1,'=,E', & I2,'.',I2,')') IF(IUNIT.EQ.0)THEN IF (I.EQ.1) WRITE (*,FORM) R1 IF (I.EQ.2) WRITE (*,FORM) R2 ELSE IF (I.EQ.1) WRITE (IUNIT,FORM) R1 IF (I.EQ.2) WRITE (IUNIT,FORM) R2 ENDIF 24 CONTINUE IF (LKNTRL.LE.0) GO TO 40 ! ERROR NUMBER IF(IUINT.EQ.0)THEN WRITE(*,30) LERR ELSE WRITE (IUNIT,30) LERR ENDIF 30 FORMAT (15H ERROR NUMBER =,I10) 40 CONTINUE 50 CONTINUE ! TRACE-BACK IF (LKNTRL.GT.0) CALL FDUMP 100 CONTINUE IFATAL = 0 IF ((LLEVEL.EQ.2).OR.((LLEVEL.EQ.1).AND.(MKNTRL.EQ.2))) & IFATAL = 1 ! QUIT HERE IF MESSAGE IS NOT FATAL IF (IFATAL.LE.0) RETURN IF ((LKNTRL.LE.0).OR.(KOUNT.GT.MAX0(1,MAXMES))) GO TO 120 ! PRINT REASON FOR ABORT IF (LLEVEL.EQ.1) CALL XERPRT & ('JOB ABORT DUE TO UNRECOVERED ERROR.',35) IF (LLEVEL.EQ.2) CALL XERPRT & ('JOB ABORT DUE TO FATAL ERROR.',29) ! PRINT ERROR SUMMARY CALL XERSAV(' ',-1,0,0,KDUMMY) 120 CONTINUE ! ABORT IF ((LLEVEL.EQ.2).AND.(KOUNT.GT.MAX0(1,MAXMES))) LMESSG = 0 CALL XERABT(MESSG,LMESSG) RETURN END SUBROUTINE XERSAV(MESSG,NMESSG,NERR,LEVEL,ICOUNT) !***BEGIN PROLOGUE XERSAV !***DATE WRITTEN 800319 (YYMMDD) !***REVISION DATE 820801 (YYMMDD) !***CATEGORY NO. Z !***KEYWORDS ERROR,XERROR PACKAGE !***AUTHOR JONES, R. E., (SNLA) !***PURPOSE Records that an error occurred. !***DESCRIPTION ! From the book "Numerical Methods and Software" ! by D. Kahaner, C. Moler, S. Nas