gistfile1.f90

PROGRAM example
REAL enought,t,phi
	
WRITE(6,*)'Enter accelerating voltage in kV'
READ(5,*)enought
WRITE(6,*)'Enter sample thickness in Angstroms'
READ(5,*)t
WRITE(6,*)'Enter cone semi-angle in mrad'
READ(5,*)phi
CALL prec_2beam(enought,t,phi)
WRITE(6,*) 'Done.'
END   

gistfile2.f90

c Subroutine to do a forward-calculation of 2-beam precession intensity.

c This is an integration around the precession circuit of the equation 
c for 2-beam intensity (see e.g. Hirsch p. 202).  

c Variable definitions
c Input: enought (accelerating voltage in kV), t (sample thickness in Angstroms),phi (cone semi-angle in mrad)
C Output: '2beam_intensity.txt'
c xlam  - electron wavelength in vacuum, in Angstroms (adjusted for mean inner potential) 
c t     - sample thickness in Angstroms 
c phi   - precession cone semi-angle in radians
c u     - structure factor U(g) for the reflection of interest (g) (in A^-2) 
c g     - spatial frequency (1/d) for reflection g (in A^-1)
c pi    - value of pi=3.14159265
c------------------------------------------------------------------------------------
c Subroutine computes the 2-beam precession intensity, deriv wrt to u, deriv wrt t for each hkl to output file

	SUBROUTINE prec_2beam(enought,t,phi)
C	implicit double precision(A-Z)
	
	REAL enought,xlam,t,phi,u,g,pi
	REAL xint,xintA,xintB,dIdu,dIdt,ustep,tstep,u0,t0
	EXTERNAL func     
	PARAMETER (maxref=1000)
	INTEGER ih(maxref),k(maxref),l(maxref)
	REAL uc(maxref),gc(maxref)
	
	pi=3.141592654
	enought=enought*1000.
	xk=sqrt(enought*(1.0+.97845E-6*enought))/12.2639
	xlam=1./xk
	phi=phi/1000.
	
c read in reflections for GITO to get h,k,l, |u| and |g|: Input file should list
c the fourier coefficient of the potential, |V(g)| which is independent of the
c accelerating voltage.  

c value of 2m/h**2:
	uconv=6.648352E-3*(1.E0+1.956934E-6*enought)

	open(1,file='gito_spots.txt')
	read(1,*)nspots

	DO i=1,nspots
		read(1,*)ih(i),k(i),l(i),vc,gc(i)
c convert here from V(g) to U(g):
		uc(i)=vc*uconv
		if(ih(i).eq.0.and.l(i).eq.0)unought=uc(i)
	ENDDO
	close(1)
	
c correct the wavelength for mean inner potential:
	xk=xk+unought/(2.*xk)
	xlam=1./xk
	
	open(2,file='2beam_intensity.txt')

	DO i=1,nspots
		IF(ih(i).eq.0.and.k(i).eq.0.and.l(i).eq.0)THEN
c zero beam not valid - skip it.  
			xint=0.
			dIdu=0.
			dIdt=0.
		ELSE
			g=gc(i)
			u=uc(i)
			xint=0.
			xintA=0.
			xintB=0.
  	
  	
c evaluate the integral around the precession circuit (over angle alpha) using qtrap.  

c func is a function which provides the integrand.  Limits are taken as 0 ... pi
c because by the symmetry of the problem we can do this then multiply by 2.  

			call qtrap(func,0.0,pi,xint,phi,u,g,pi,xlam,t)

c multiply by 2*prefactor; divide by 2*pi (effectively integrating over 2*pi): 
			xint=2.*xint*(pi*xlam*u)**2/(2.*pi)

			u0=u
			t0=t
			ustep=1d-4
			tstep=1d-4 !accurate to 8 decimal places
c Numerical derivative dI/du	
			u=ustep+u0
			call qtrap(func,0.0,pi,xintA,phi,u,g,pi,xlam,t)
			xintA=2.*xintA*(pi*xlam*u)**2/(2.*pi)
			u=u0-ustep
			call qtrap(func,0.0,pi,xintB,phi,u,g,pi,xlam,t)
			xintB=2.*xintB*(pi*xlam*u)**2/(2.*pi)
			dIdu=(xintA-xintB)/(2.*ustep)
C Numerical derivative dI/dt
			u=u0
			t=t0+tstep
			call qtrap(func,0.0,pi,xintA,phi,u,g,pi,xlam,t)
			xintA=2.*xintA*(pi*xlam*u)**2/(2.*pi)
			t=t0-tstep
			call qtrap(func,0.0,pi,xintB,phi,u,g,pi,xlam,t)
			xintB=2.*xintB*(pi*xlam*u)**2/(2.*pi)
			dIdt=(xintA-xintB)/(2.*tstep)
			t=t0
		ENDIF
C Write to output file '2beam_intensity.txt'
		write(2,'(3i5,3e17.10)')ih(i),k(i),l(i),xint,dIdu,dIdt
	ENDDO   
	                                    
	close(2)
     
	RETURN 
	END    
         
c**************************************************************
c 2-beam sinc function, with s(eff) as in Hirsch eq. 8.23.  
c Evaluate here without the leading term.   
c alpha is the precession azimuth angle (over which we're integrating)
	
	REAL FUNCTION func(alpha,phi,u,g,pi,xlam,t)
	REAL xlam,t,phi,u,g,pi,alpha,arg,s
	
c get s as function of g, phi and alpha:  
	s=-xlam*g**2/2.+g*phi*cos(alpha)

c s(effective) for 2-beam.  For kinematical calculation, set sef=s.  
	sef=sqrt(s*s+xlam*xlam*u*u)
c	sef=s
	
	arg=pi*t*sef
	IF(abs(arg).lt.1.e-5)THEN
		func=t*t
	ELSE
		func=sin(arg)/(pi*sef)
		func=func*func
	ENDIF
	
	RETURN
	END

c ********* NUMERICAL RECIPES CODES ***************
c*******************************************************
      SUBROUTINE qtrap(func,a,b,s,phi,u,g,pi,xlam,t)
      INTEGER JMAX,j
      REAL a,b,func,s,EPS,phi,u,g,pi,xlam,t,olds
      EXTERNAL func
      PARAMETER (EPS=2.e-6, JMAX=20)
C uses trapzd

      olds=-1.e30
      DO 11 j=1,JMAX
        call trapzd(func,a,b,s,j,phi,u,g,pi,xlam,t)
        IF (abs(s-olds).lt.EPS*abs(olds)) RETURN
        IF (s.eq.0..and.olds.eq.0..and.j.gt.6) RETURN
        olds=s
11    CONTINUE
      write(6,*)'too many steps in qtrap'
      RETURN
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 510.
c*****************************************************************
      SUBROUTINE trapzd(func,a,b,s,n,phi,u,g,pi,xlam,t)
      INTEGER n,it,j
      REAL a,b,func,s,phi,u,g,pi,xlam,t,del,sum,tnm,x
      EXTERNAL func
      
      IF (n.eq.1) THEN
       s=0.5*(b-a)*(func(a,phi,u,g,pi,xlam,t)+func(b,phi,u,g,pi,xlam,t))
      ELSE
        it=2**(n-2)
        tnm=it
        del=(b-a)/tnm
        x=a+0.5*del
        sum=0.
        DO 11 j=1,it
          sum=sum+func(x,phi,u,g,pi,xlam,t)
          x=x+del
11      CONTINUE
        s=0.5*(s+(b-a)*sum/tnm)
      ENDIF
      RETURN
      END