!  Superflourescence in laser amplifiers
!  Source file:  laser.f90		header:  laser.h
!
!					Ernest Sibert - November 1996

!  Declared or defined in laser.h:
!
!	mu	index of refraction
!	alpha	amplification coefficient
!	psicrit	critical angle (radians) = asin(1/mu)
!	r	radius of disk
!	t	thickness of disk
!	ht	t/2
!	pi	3.14159...

!  Array organization:
!
!	real, array(1:Nx, 0:Nz, 0:Nphi) corresponding to x, z, phi

!  Calculate the integral using specified array dimensions, steps
!  for psi-regions.

      real function superf(Nx, Nz, Nphi, N1, N2, N3)
      integer Nx, Nz, Nphi, N1, N2, N3

      INCLUDE 'laser.h'

      real, dimension(1:Nx, 0:Nz, 0:Nphi) ::             &
	x, z, phi, psi, lng, h, v, f, integ, w1, w2, w4, &
	hpsi, psi1, psi2
      integer, dimension(1:Nx, 0:Nz, 0:Nphi) :: Ix, Iz, Iphi
      real  hx, hz, hphi
      integer ips, NN

      NN   = Nx*Nz*Nphi		! Array size
      hx   = r/Nx		! Interval in x
      hz   = ht/Nz		! Interval in z
      hphi = pi/Nphi		! Interval in phi
!     print *, "hx, hz, hphi: ", hx, hz, hphi
      Ix   = spread(spread([1:Nx],   2, Nz+1), 3, Nphi+1)	! Coordinates
      Iz   = spread(spread([0:Nz],   1, Nx),   3, Nphi+1)
      Iphi = spread(spread([0:Nphi], 1, Nx),   2,  Nz+1)
      x    = hx*Ix
      where ( Ix == Nx ) x = r		! On account of rounding errors
      z    = hz*Iz
      where ( Iz == Nz ) z = ht
      phi  = hphi*Iphi
      where ( Iphi == Nphi ) phi = pi
      h    = sqrt(r**2 - (x*sin(phi))**2) - x*cos(phi)	! Distance to side
      where ( h <= 0.0 ) h = 0.0
!     print *, "h, min, avg, max ", &
!    		minval(h), sum(h)/NN, maxval(h)
      v    = ht - z					! Distance to top
      where ( z < ht )		! All except top
	psi1 = min( psicrit, atan(h/v) )
      elsewhere
	psi1 = 0.5*pi
      end where
!     print *, "psi1, min, avg, max ", &
!    		minval(psi1), sum(psi1)/NN, maxval(psi1)

!  psi integral from 0 to psi1

      hpsi = psi1/N1
!     print *, "hpsi, min, avg, max ", &
!    		minval(hpsi), sum(hpsi)/NN, maxval(hpsi)
      w2 = 0.0
      w4 = 0.0
      do, ips = 1, N1-1				! Integrand is 0 for psi = 0
	psi = ips*hpsi
	f = exp(alpha*v/cos(psi))*sin(psi)	! Integrand
        if ( iand(ips, 1) == 1 ) then		! ips is odd
	  w4 = w4 + f
	else					! ips is even
	  w2 = w2 + f
	end if
      end do ! iph = 1, N1-1
      f  = 1.0					! For z = t/2
      where ( Iz < Nz)	&			! Ill-defined on top
        f = exp(alpha*v/cos(psi1))*sin(psi1)	! last integrand
      integ = hpsi*(f + 4.0*w4 + 2.0*w2)/3.0
!     print *, "int1, min, avg, max", &
!    		minval(integ), sum(integ)/NN, maxval(integ)

!  psi integral from psi1 to psi2

      v = ht + z				! Distance to bottom
      psi2 = max( pi - psicrit, pi - atan(h/v) )
!     print *, "psi2, min, avg, max ", &
!    		minval(psi2), sum(psi2)/NN, maxval(psi2)
      hpsi = (psi2 - psi1)/N2
!     print *, "hpsi, min, avg, max ", &
!    		minval(hpsi), sum(hpsi)/NN, maxval(hpsi)
      w1 = 0.0
      w2 = 0.0
      w4 = 0.0
      f = sin(psi1)
      where ( h > 0.0 ) &			! Ill-formed at sidewall
        w1 = exp(alpha*h/f)*f			! Integrand at psi1
      where ( h == 0.0 .and. Iz == Nz ) &	! Top and sidewall
      	w1 = 1.0				! psi1 = pi/2, lng = 0
      f = sin(psi2)
      where ( h > 0.0 ) &			! Ill-formed at sidewall
        w1 = w1 + exp(alpha*h/f)*f		! Integrand at psi2
      do, ips = 1, N2-1				! Intermediate points
	psi = psi1 + ips*hpsi
	f = sin(psi)
	f = exp(alpha*h/f)*f
        if ( iand(ips, 1) == 1 ) then		! ips is odd
	  w4 = w4 + f
	else					! ips is even
	  w2 = w2 + f
	end if
      end do ! ips = 1, N2-1
      integ = integ + hpsi*(w1 + 4.0*w4 + 2.0*w2)/3.0
!     print *, "+ int2, min, avg, max", &
!    		minval(integ), sum(integ)/NN, maxval(integ)
!  psi integral from psi2 to pi
      hpsi = (pi - psi2)/N3
!     print *, "hpsi, min, avg, max ", &
!    		minval(hpsi), sum(hpsi)/NN, maxval(hpsi)
      w1 = exp(-alpha*v/cos(psi2))*sin(psi2)
      w2 = 0.0
      w4 = 0.0
      do, ips = 1, N3-1				! Integrand is 0 at psi = pi
	psi = psi2 + ips*hpsi
	f = exp(-alpha*v/cos(psi))*sin(psi)
        if ( iand(ips, 1) == 1 ) then		! ips is odd
	  w4 = w4 + f
	else					! ips is even
	  w2 = w2 + f
	end if
      end do ! ips = 1, N3-1
      integ = integ + hpsi*(w1 + 4.0*w4 + 2.0*w2)/3.0
!     print *, "+ int3, min, avg, max", &
!    		minval(integ), sum(integ)/NN, maxval(integ)
!  Now integ = complete psi integral
      f = x*integ				! Outer integrand
!     print *, "f, min, avg, max ", &
!    		minval(f), sum(f)/NN, maxval(f)
!  Compute weights for 3-D Simpson
      w1 = 1.0
      where ( iand(Ix,1) == 1 )  w1 = 4.0
      where ( iand(Ix,1) == 0 .and. Ix < Nx ) w1 = 2.0
      where ( iand(Iz,1) == 1 ) w1 = 4.0*w1
      where ( iand(Iz,1) == 0 .and. 0 < Iz .and. Iz < Nz) &
	w1 = 2.0*w1
      where ( iand(Iphi,1) == 1) w1 = 4.0*w1
      where ( iand(Iphi,1) == 0 .and. 0 < Iphi .and. Iphi < Nphi ) &
	w1 = 2.0*w1
!     print *, "w1, min, avg, max ", &
!    		minval(w1), sum(w1)/NN, maxval(w1)
      superf = 2.0*hx*hz*hphi*SUM(w1*f)/27.0
      end function superf


!  Program to manage calculations

      program laseramp
      INCLUDE 'laser.h'

      INTERFACE
      real function superf(Nx, Nz, Nphi, N1, N2, N3)
      integer Nx, Nz, Nphi, N1, N2, N3
      END INTERFACE

      real PHIovN
      integer Nx, Nz, Nphi, N1, N2, N3

      print *, "Laser material parameters:"
      print '(" mu alpha > "$)'
      read *, mu, alpha
      psicrit = asin(1.0/mu)
      print 21, psicrit, psicrit*180.0/pi
 21   FORMAT (" Critical angle = ", F8.6, " (", F8.5, " deg.)")
!  Loop for laser dimension
      do while ( .true. )
 99   print *, "Disk dimensions ( r <= 0 to stop ):"
      print '(" r t > "$)'
      read *, r, t
      if ( r <= 0.0 ) stop			! Stop
      ht = 0.5*t
      print *, "Integration parameters, all must be even."
      print *, "Specify Nx <= 0 to exit loop."
      do while ( .true. )
 88	print '(" Nx Nz Nphi > "$)'
        read *, Nx, Nz, Nphi
	if ( Nx <= 0 ) go to 99		! New dimensions
	print '(" N1  N2  N3 > "$)'
	read *, N1, N2, N3
	if ( &		! Any parameter not even
	    iand(ior(Nx, ior(Nz, ior(Nphi, ior(N1, ior(N2,N3))))),1) > 0 &
	   ) then
	  print *, "Integration parameters MUST be EVEN."
	  go to 88
	end if
	PHIovN = superf(Nx, Nz, Nphi, N1, N2, N3)
        print *, "PHI/N = ", PHIovN, &
		 "  PHI/(N*V) = ", PHIovN/(pi*r*r*t)
      end do ! Integration
      end do ! Dimensions
      end program laseramp
