text:
Dave can describe himself here.  He appears to be Dave Reynolds of Paintings by Dave.

20110413 14:06:07 nbsp Hi Dave, and Welcome to the Wiki!  Take a moment to look over the Welcome to the Wiki/Business Owner business owner welcome pagetheres a lot of very useful info there.  I touched up the table on your new page a little bithopefully it looks OK.  Thanks for adding your business to the wiki!  You can write up a blurb about Paintings by Dave there, and maybe include some info about the artist (yourself).  One great thing about the wiki is learning a bit about the local folks we meet and/or do business with.  See that business owner welcome page for tips on how it should lookit should be informational in nature, so it reads as being information about you and the business, rather than from you and the business.  The next step would be to link your business from some related pages, and to add links to related stuff on the Paintings by Dave page.  For example, if its an art studio, you might want to add a link (preferably in alphabetical order) to the Art page.  Its quite easy to do  just hit the edit button at the top of that page, find the right area, and add in {{{  Paintings by Dave}}}.  Other pages that might be relevant include Art Galleries, E Street, and Downtown. Users/TomGarberson

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      double precision function  fdisc(blockno,xc,yc)
      implicit none

      double precision xc,yc, rc,xloc(0:4),yloc(0:4)
      integer blockno
      integer i

      data xloc /0, 1, 1, 0, 0.5d0/
      data yloc /0, 0, 1, 1, 0.5d0/

c     # Map each brick to a [0,1]x[0,1] domain and duplicate
c     # initial conditions.
      do i = 0,4
         rc = sqrt((xc-xloc(i))**2 + (yc-yloc(i))**2)
         fdisc = rc-0.3d0
         if (fdisc .lt. 0) then
            return
         endif
      enddo

      end

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


text:
!     .................................................
!             ____  _       _   ____  _____   _
!            |  _ \| |     |_| |  _ \|  ___| |_|
!            | |_) | |___   _  | |_) | |___   _
!            |  _ /|  _  | | | |  _ /|___  | | |
!            | |   | | | | | | | |    ___| | | |
!            |_|   |_| |_| |_| |_|   |_____| |_|
!     .................................................
!     PhiPsi:     a general-purpose computational
!                 mechanics program written in Fortran.
!     Website:    http://phipsi.top
!     Author:     Shi Fang from Huaiyin Institute of
!                 Technology, HuaiAn, JiangSu, China
!     Contact me: shifang@hyit.edu.cn
!     ------------------------------------------------
!     Please cite the following papers:
!     (1)Shi F, Wang X L, Liu C, Liu H, Wu H A. An
!        XFEM-based method with reduction technique
!        for modeling hydraulic fracture propagation
!        in formations containing frictional natural
!        fractures. Engineering Fracture Mechanics,
!        2017, 173: 64-90.
!     (2)Shi F, Wang X L, Liu C, Liu H, Wu H A. A
!        coupled extended finite element approach
!        for modeling hydraulic fracturing in
!        consideration of proppant. Journal of
!        Natural Gas Science and Engineering, 2016,
!        33: 885-897.
!     (3)Shi F, Wang X L, Liu C, Liu H, Wu H A. An
!        XFEM-based numerical model to calculate
!        conductivity of propped fracture considering
!        proppant transport, embedment and crushing.
!        Journal of Petroleum Science and Engineering,
!        2018, 167: 615-626..

      SUBROUTINE Save_Disp(isub)
c     save displacement.

      use Global_Common
      use Global_Model
      use Global_Filename
      use Global_DISP

      implicit none

      integer isub
      integer i
      character(200) c_File_name_1
      character(5) temp

      print *,'    Saving disp......'
      write(temp,'(I4)') isub
      c_File_name_1   =  trim(Full_Pathname)//'.'
     &                 //'disp'//'_'//ADJUSTL(temp)

      open(101,file=c_File_name_1,status='unknown')
      do i=1,Total_Freedom
          write(101, '(1E20.12)') DISP(i)
      end do
      close(101)

      write(temp,'(I4)') isub
      c_File_name_1   =  trim(Full_Pathname)//'.'
     &             //'disn'//'_'//ADJUSTL(temp)
      open(101,file=c_File_name_1,status='unknown')
      do i=1,Total_Freedom
          write(101, '(I8,2E20.12)') i,DISP(2*i-1),DISP(2*i)
      end do
      close(101)

      RETURN
      END SUBROUTINE Save_Disp


text:
! This program calculates the matrix elements of the positronium-hydrogen
!  matrices (similar to equation 2.15 of the Armour and Humberston article).  Specifically,
!  we are calculating elements of the form (phi_i, L phi_j) as in equation (3.22).
! This does use OpenMP to speed up computation on multicore processors, so some familiarity
!  with parallel programming is required to understand the !$omp lines.
! A future speed-up could be done with Open MPI (as is done for the energy code), but this
!  computation is usually done in a few hours.

! TODO: Make MaxIter and the 4 come from the configuration file.

module WLimits
   implicit none
   integer lmin, lmax
   integer mmin, mmax
   integer nmin, nmax
end module WLimits

program PsHMain
	use WLimits
	implicit none

	interface
		integer(c_int) function optimizewavefn(MaxIter, Omega, NumTerms, M12max, M23max, M31max, pmax, IsTriplet, Ordering, Method, EigenRoutine, alpha, beta, gamma) bind(c, name='optimizewavefn')
		use iso_c_binding
		integer(c_int), value :: MaxIter, Omega, NumTerms, M12max, M23max, M31max, pmax, IsTriplet, Ordering, Method, EigenRoutine
		real(c_double), value :: alpha, beta, gamma
		end function optimizewavefn
	end interface

	real*16 omp_get_wtime
	integer Omega  ! This sets the limits on the terms.
	integer NumTerms  ! Number of terms in our Hylleraas expansion
	integer iread, iwrite  ! Our input and output files
	integer CalcPowerTableSize
	real*16, allocatable, dimension(:,:) :: PhiPhi, Phi2HPhi
	real*8, allocatable, dimension(:,:) :: PhiPhi8, Phi2HPhi8
	real*8, allocatable, dimension(:) :: Energies, Workspace
	real*16 Alpha, Beta, Gamma, Prev
	integer, allocatable, dimension(:) :: UsedTerms
	integer NumUsed, i, j, k, n
	real*8 StartTime, EndTime
	integer iargc, IsTriplet, Method, LowerEigen, UpperEigen, Info, AllEnergies, EigenRoutine
	character *100 IOBuffer
	integer Optimize, EigenNum, Ordering, Iter, MaxIter, MaxIterOuter, M12max, M23max, M31max, pmax, StartEnergy, LValue
	real*16 Tol, Err, h, EDeriv, fa, fh, y, Divisor, ederivwrapper
	real*16, dimension(3) :: Del, Delta, f1, f0, x1, x0, x, C31, C13
	real*16, dimension(3,3) :: B, BTemp, C33, Jacob
    real*16 ThreeJSymbol, HylleraasIntegralGeneral

!	real*16, allocatable, dimension(:,:,:,:) :: WMatrix
!	call GenFactorial(400)
!	lmin = 1;  lmax = 1
!	mmin = 1;  mmax = 1
!	nmin = 1;  nmax = 1
!	allocate(WMatrix(lmin:lmax, mmin:mmax, nmin:nmax, 6))
!    write (*,*) HylleraasIntegralGeneral(.true., .false., 2, 2, 0, -1, 1, 1, 1.0q0, 1.2q0, 2.2q0, &
!								2, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 10, 150, 5, 51, 0, 0, 0, WMatrix, 1)
!    deallocate(WMatrix)

	iread = 9
	iwrite = 10
	! This allows the possibility of using different input and output files than the defaults.
	!  This mainly lets me run this program on multiple computers at the same time without
	!  one overwriting the results from another.  If the program is called without any
	!  command-line arguments, it uses the default filenames.
	if (iargc() == 2) then
		call getarg(1, IOBuffer)
		open(iread, FILE=IOBuffer)
		call getarg(2, IOBuffer)
		open(iwrite, FILE=IOBuffer)
	else
		open(iread, FILE='input.txt')
		open(iwrite, FILE='output.txt')
	endif

	! Read the input parameters.
	call ReadParamFile(iread, Omega, LValue, Alpha, Beta, Gamma, M12max, M23max, M31max, pmax, Method, IsTriplet, LowerEigen, UpperEigen, &
						 AllEnergies, Optimize, EigenNum, MaxIter, Ordering, EigenRoutine)

	if (IsTriplet /= 0 .and. IsTriplet /= 1) then
		write (*,*) 'Must choose between singlet and triplet calculation - exiting.'
		stop
	endif

	if (Optimize == 1) write (iwrite,*) 'Optimizing eigenvalue', EigenNum

	Tol = 1e-3
	Err = 1  ! Just has to be larger than tol.
	h = 1e-5

	! Get the start time to determine the duration of the program.
	!call cpu_time(StartTime)
	StartTime = omp_get_wtime()

	NumTerms = CalcPowerTableSize(Omega) * 2  ! First and second symmetries
	call WriteHeader(iwrite, LValue, IsTriplet, Ordering, Method, EigenRoutine, Omega, NumTerms, Alpha, Beta, Gamma, M12max, M23max, M31max, pmax)

	write (*,*) 'NumTerms:', NumTerms
	allocate(PhiPhi(NumTerms,NumTerms))
	allocate(Phi2HPhi(NumTerms,NumTerms))
	allocate(PhiPhi8(NumTerms,NumTerms))
	allocate(Phi2HPhi8(NumTerms,NumTerms))
	allocate(Energies(NumTerms))
	allocate(Workspace(3*NumTerms-1))

	if (Omega == -1) then
		goto 200  ! No short-range terms, so no need to do these calculations (just a header)
	end if

	if (Optimize == 1) then
		! TODO: Add these back in.
		!fa = optimizewavefn(MaxIter, Omega, NumTerms, M12max, M23max, M31max, pmax, IsTriplet, Ordering, Method, EigenRoutine, Alpha, Beta, Gamma)
		!fa = EDeriv(Omega, NumTerms, Alpha, Beta, Gamma, EigenNum, PhiPhi, Phi2HPhi, M12max, pmax, 0, &
		!			IsTriplet, Ordering, Method)
	else
		! If we are not optimizing eigenvalues, we can still call the same function, which exits early with PhiPhi and Phi2HPhi filled.
		fa = EDeriv(Omega, NumTerms, Alpha, Beta, Gamma, EigenNum, PhiPhi, Phi2HPhi, M12max, M23max, M31max, pmax, 0, &
					IsTriplet, Ordering, Method, EigenRoutine, LValue)
		write (*,*) fa
    endif

    if (LValue == 0) then  ! S-wave only has one symmetry
        NumTerms = NumTerms / 2
    end if

	do i = 1, NumTerms, 1
		!do j = 1, i, 1  ! Use this instead if only the lower triangle is required.
		do j = 1, NumTerms, 1
			write (iwrite,"(i6,i6,d38.30,d38.30)") i, j, PhiPhi(i,j), Phi2HPhi(i,j)
		enddo
	enddo

	! Divide every entry in Phi2HPhi by 2, since we calculated <phi_i|2H|phi_j> in equation (3.22).
	Phi2HPhi = Phi2HPhi / 2.0_16

	write (iwrite,*)
	write (iwrite,*) 'Energies:'

	if (AllEnergies == 0) then
        StartEnergy = NumTerms
	else
		StartEnergy = 1
    end if

    if (LValue == 0) then
        NumTerms = NumTerms * 2
    end if

	do j = StartEnergy, NumTerms, 1
	    if (LValue == 0 .and. j > NumTerms / 2) then  ! S-wave only has one symmetry
            exit
        end if
		! Copies our real*16 matrices to real*8 matrices so that LAPACK can operate on them.
		Phi2HPhi8 = Phi2HPhi
		PhiPhi8 = PhiPhi

		! This calculates the energy eigenvalues of the generalized eigenvalue problem generated from
		!  the Rayleigh-Ritz variational method:
		!  det(<phi_i|H|phi_j> - E <phi_i|phi_j>) = 0
		! We calculated the lower and upper triangles of the PhiPhi and Phi2HPhi matrices above, though
		!  they are symmetric.  LAPACK does not require us to fill in both halves; we can just
		!  specify 'L' in the third parameter to denote that the lower triangle (and diagonal) is filled.
		! Explanation of choices for the parameters are found at the Netlib site:
		!  http://www.netlib.org/lapack/double/dsygv.f
		call dsygv(1, 'N', 'L', j, Phi2HPhi8, NumTerms, PhiPhi8, NumTerms, Energies, Workspace, 3*NumTerms-1, Info)
		!call Newdsygv(1, 'N', 'L', j, Phi2HPhi, NumTerms, PhiPhi, NumTerms, Energies, Workspace, 3*NumTerms-1, Info)
		!call nag_sym_gen_eig_all('L', Phi2HPhi8, PhiPhi8, Energies)  ! NAG equivalent
		if (Info /= 0) then
			write (*,*)
			write (*,*) 'Energy eigenvalues could not be determined!'
			write (*,*) 'dsygv error code:', Info
			write (iwrite,*)
			write (iwrite,*) 'Energy eigenvalues could not be determined!'
			write (iwrite,*) 'dsygv error code:', Info
			exit  ! Exit the loop, because every successive j will give the same error, since we have a troublesome term.
		elseif (Info == 0) then  ! dsygv successfully completed.
			! Writes the eigenvalues.
			! This output is formatted properly for inclusion into Excel as a comma-separated value (.csv) file.
			write (iwrite,"(i4)",advance='no') j
			do i = LowerEigen, UpperEigen, 1
				write (iwrite,"(a,d21.14)",advance='no') ', ', Energies(i)
			enddo
			write (iwrite,*) ' '  ! Finish the line

			write (*,"(i4)",advance='no') j
			do i = LowerEigen, UpperEigen, 1
				write (*,"(a,d21.14)",advance='no') ', ', Energies(i)
			enddo
			write (*,*) ' '  ! Finish the line
		endif
	enddo

	! Get the end time to find the duration of the program.
	!call cpu_time(EndTime)
200	EndTime = omp_get_wtime()
	write (*,*) 'Time taken (s):', EndTime - StartTime
	write (iwrite,*)
	write (iwrite,*) 'Time taken (s):', EndTime - StartTime, '  (min):', (EndTime - StartTime) / 60.0

	! Clean up memory before exiting
	deallocate(Workspace)
	deallocate(Energies)
	deallocate(PhiPhi8)
	deallocate(Phi2HPhi8)
	deallocate(PhiPhi)
	deallocate(Phi2HPhi)

	! Close file handles
	close(iwrite)

	stop
end

subroutine ReadParamFile(iread, Omega, LValue, Alpha, Beta, Gamma, M12max, M23max, M31max, pmax, Method, IsTriplet, LowerEigen, UpperEigen, &
						 AllEnergies, Optimize, EigenNum, MaxIter, Ordering, EigenRoutine)
	implicit none
	integer iread, LValue, Omega, Method, M12max, M23max, M31max, pmax, IsTriplet, LowerEigen, UpperEigen, AllEnergies
	integer Optimize, EigenNum, MaxIter, Ordering, EigenRoutine
	real*16 Alpha, Beta, Gamma

	read (iread,*) ! Description line
	read (iread,*) Omega
	read (iread,*) ! Description line
	read (iread,*) LValue
	read (iread,*)
	read (iread,*) Alpha
	read (iread,*) Beta
	read (iread,*) Gamma
	read (iread,*)
	read (iread,*) Method
	read (iread,*)
	read (iread,*) M12max, M23max, M31max
	read (iread,*)
	read (iread,*) pmax
	read (iread,*)
	read (iread,*) IsTriplet
	read (iread,*)
	read (iread,*) LowerEigen, UpperEigen
	read (iread,*)
	read (iread,*) AllEnergies
	read (iread,*)
	read (iread,*) Optimize
	read (iread,*)
	read (iread,*) EigenNum
	read (iread,*)
	read (iread,*) MaxIter
	read (iread,*)
	read (iread,*) Ordering
	read (iread,*)
	read (iread,*) EigenRoutine
	close (iread)

	return
end

subroutine WriteHeader(iwrite, LValue, IsTriplet, Ordering, Method, EigenRoutine, Omega, NumTerms, Alpha, Beta, Gamma, M12max, M23max, M31max, pmax)
	implicit none
	integer iwrite, LValue, IsTriplet, Ordering, Method, EigenRoutine, Omega, NumTerms, M12max, M23max, M31max, pmax
	real*16 Alpha, Beta, Gamma

	select case (LValue)
		case (0)  ! S-Wave
			if (IsTriplet == 0) then
				write (*,*) "S-Wave Singlet Ps-H"
				write (iwrite,*) "S-Wave Singlet Ps-H"
			else
				write (*,*) "S-Wave Triplet Ps-H"
				write (iwrite,*) "S-Wave Triplet Ps-H"
			endif
		case (1)
			if (IsTriplet == 0) then
				write (*,*) "P-Wave Singlet Ps-H: 1st formalism"
				write (iwrite,*) "P-Wave Singlet Ps-H: 1st formalism"
			else
				write (*,*) "P-Wave Triplet Ps-H: 1st formalism"
				write (iwrite,*) "P-Wave Triplet Ps-H: 1st formalism"
			endif
		case (2)
			if (IsTriplet == 0) then
				write (*,*) "D-Wave Singlet Ps-H: 1st formalism"
				write (iwrite,*) "D-Wave Singlet Ps-H: 1st formalism"
			else
				write (*,*) "D-Wave Triplet Ps-H: 1st formalism"
				write (iwrite,*) "D-Wave Triplet Ps-H: 1st formalism"
			endif
		case (3)  ! F-Wave
			if (IsTriplet == 0) then
				write (*,*) "F-Wave Singlet Ps-H"
				write (iwrite,*) "F-Wave Singlet Ps-H"
			else
				write (*,*) "F-Wave Triplet Ps-H"
				write (iwrite,*) "F-Wave Triplet Ps-H"
			endif
		case (4)  ! G-Wave
			if (IsTriplet == 0) then
				write (*,*) "G-Wave Singlet Ps-H"
				write (iwrite,*) "G-Wave Singlet Ps-H"
			else
				write (*,*) "G-Wave Triplet Ps-H"
				write (iwrite,*) "G-Wave Triplet Ps-H"
			endif
		case (5)  ! H-Wave
			if (IsTriplet == 0) then
				write (*,*) "H-Wave Singlet Ps-H"
				write (iwrite,*) "H-Wave Singlet Ps-H"
			else
				write (*,*) "H-Wave Triplet Ps-H"
				write (iwrite,*) "H-Wave Triplet Ps-H"
			endif
		case (6)  ! I-Wave
			if (IsTriplet == 0) then
				write (*,*) "I-Wave Singlet Ps-H"
				write (iwrite,*) "I-Wave Singlet Ps-H"
			else
				write (*,*) "I-Wave Triplet Ps-H"
				write (iwrite,*) "I-Wave Triplet Ps-H"
			endif
		case (7)  ! K-Wave
			if (IsTriplet == 0) then
				write (*,*) "K-Wave Singlet Ps-H"
				write (iwrite,*) "K-Wave Singlet Ps-H"
			else
				write (*,*) "K-Wave Triplet Ps-H"
				write (iwrite,*) "K-Wave Triplet Ps-H"
			endif
		case (8)  ! L-Wave
			if (IsTriplet == 0) then
				write (*,*) "L-Wave Singlet Ps-H"
				write (iwrite,*) "L-Wave Singlet Ps-H"
			else
				write (*,*) "L-Wave Triplet Ps-H"
				write (iwrite,*) "L-Wave Triplet Ps-H"
			endif
		case default
			write (*,*) "Higher partial waves are not supported yet...exiting."
			stop
	end select

	if (Ordering == 1) then
		write (*,*) "Using Peter Van Reeth's ordering"
		write (iwrite,*) "Using Peter Van Reeth's ordering"
	else
		write (*,*) "Using my ordering"
		write (iwrite,*) "Using my ordering"
	endif

	if (Method == 0) then
		write (*,*) "Integration Technique: Direct Summation"
		write (iwrite,*) "Integration Technique: Direct Summation"
	else if (Method == 1) then
		write (*,*) "Integration Technique: Asymptotic Expansion"
		write (iwrite,*) "Integration Technique: Asymptotic Expansion"
	else if (Method == 2) then
		write (*,*) "Integration Technique: Recursion Relations"
		write (iwrite,*) "Integration Technique: Recursion Relations"
	else
		write (*,*) "Method parameter in input file must be 0, 1 or 2."
		stop
	end if

	if (EigenRoutine == 1) then
		write (*,*) "Eigenvalue routine: dsygv from LAPACK"
		write (iwrite,*) "Eigenvalue routine: dsygv from LAPACK"
	else if (EigenRoutine == 2) then
		write (*,*) "Eigenvalue routine: Pachucki eigenproblem solver"
		write (iwrite,*) "Eigenvalue routine: Pachucki eigenproblem solver"
	else
		write (*,*) "Eigenvalue routine parameter in input file must be 1 or 2."
		stop
	end if

	write (iwrite,*) 'Omega =', Omega
	write (iwrite,*) 'Alpha =', Alpha
	write (iwrite,*) 'Beta =', Beta
	write (iwrite,*) 'Gamma =', Gamma

    if (LValue == 0) NumTerms = NumTerms / 2
	write (iwrite,*) 'Number of terms:', NumTerms
	write (*,*) 'Number of terms:', NumTerms
	write (*,*) 'Omega =', Omega
    if (LValue == 0) NumTerms = NumTerms * 2

	write (iwrite,*) 'M12max, M23max, M31max, pmax =', M12max, M23max, M31max, pmax
	write (iwrite,*)

	return
end

! This subroutine calculates the number of terms for the function GenOmegaPowerTable (below).  This
!  is needed to determine the size of the dynamically allocated table before it is produced.
integer function CalcPowerTableSize(Omega)
	implicit none
	integer Omega  ! This sets the limits on the terms.
	integer NumTerms  ! The total number of terms
	integer om, ki, li, mi, ni, pi, qi  ! These are the exponents we are determining.

	if (Omega == -1) then  ! Special case where we don't want any short-range terms
		CalcPowerTableSize = 0
		return
	end if

	NumTerms = 0
	do om = 0, Omega, 1
		do ki = 0, Omega, 1
			do li = 0, Omega, 1
				do mi = 0, Omega, 1
					do ni = 0, Omega, 1
						do pi = 0, Omega, 1
							do qi = 0, Omega, 1
								if (ki + li + mi + ni + pi + qi == om) then
									!if (ki>li) cycle
									!if (ki==li .and. qi>pi) cycle

!									if (om > 3 .and. qi > 0) cycle  ! Restrict r23 to max of 1
!									if (om > 3 .and. pi > 0) cycle  ! Restrict r13 to max of 3
!									if (om > 2 .and. ni > 0) cycle  ! Restrict r3 to max of 0
!									if (ki >= li .and. qi > pi) cycle  ! Check if r1>r2 and r23 > r13 (powers)
									NumTerms = NumTerms + 1
								endif
							enddo
						enddo
					enddo
				enddo
			enddo
		enddo
	enddo

	CalcPowerTableSize = NumTerms
	return
end

! This subroutine calculates values of k_i, l_i, m_i, n_i, p_i and q_i of
!  (3.14).  The summation of these is given by equation (3.15), but we do not have
!  the restriction on q of being even.
! This is written to always have the terms in order of increasing omega, i.e. the terms
!  for omega = 0 are first, etc.
! TODO: Use logical for Ordering
subroutine GenOmegaPowerTable(Omega, PowerTable, ArraySize, Ordering, l1, l2)
	!implicit none
	integer Omega  ! This sets the limits on the terms.
	integer ArraySize
	integer, dimension(ArraySize,6) :: PowerTable
	integer NumTerm  ! The number of the current term
	integer om, ki, li, mi, ni, pi, qi  ! These are the exponents we are determining.
	integer Ordering  ! Whether to use Peter's ordering or mine
	integer l1, l2, l3  ! l1 = 0 for S-wave, 1 for P-wave, etc.

	if (Ordering == 0) then  ! Use my ordering
		NumTerm = 0
		do om = 0, Omega, 1
			do ki = 0, Omega, 1
				do li = 0, Omega, 1
					do mi = 0, Omega, 1
						do ni = 0, Omega, 1
							do pi = 0, Omega, 1
								do qi = 0, Omega, 1
									if (ki + li + mi + ni + pi + qi == om) then
										!if (ki + li + mi + ni + pi + qi <= Omega) then
										!if (ki>li) cycle
										!if (ki==li .and. qi>pi) cycle

	!									if (om > 3 .and. qi > 0) cycle  ! Restrict r23 to max of 1
	!									if (om > 3 .and. pi > 0) cycle  ! Restrict r13 to max of 3
	!									if (om > 2 .and. ni > 0) cycle  ! Restrict r3 to max of 0
	!									if (ki >= li .and. qi > pi) cycle  ! Check if r1>r2 and r23 > r13 (powers)

										NumTerm = NumTerm + 1
										PowerTable(NumTerm,1) = ki
										PowerTable(NumTerm,2) = li
										PowerTable(NumTerm,3) = mi
										PowerTable(NumTerm,4) = ni
										PowerTable(NumTerm,5) = pi
										PowerTable(NumTerm,6) = qi

										write (*,*) ki, li, mi, ni, pi, qi!, ki + li + mi + ni + pi + qi
									endif
								enddo
							enddo
						enddo
					enddo
				enddo
			enddo
		enddo

	else  ! Use Peter Van Reeth's ordering instead.
		! I pulled this (modified) snippet from Peter Van Reeth's code so that
		!  my indices match up with his to debug mine easier.
		MEGA=OMEGA
			IHDPP1 = MEGA +1
			INX=0
			DO 65 I=1,IHDPP1
			  DO 66 I23P1=1,I,1
				I23=I23P1 -1
				I12P1M = I-I23
				DO 67 I12P1=1,I12P1M
				  I12=I12P1-1
				  I2P1M =I12P1M -I12
				  DO 68 I2P1 =1,I2P1M
					 I2 = I2P1-1
					 I13P1M = I2P1M-I2
					 DO 69 I13P1=1,I13P1M
						I13 = I13P1 -1
						I3P1M = I13P1M -I13
						DO 70 I3P1 = 1, I3P1M
						   I3=I3P1-1
						   I1 = I3P1M -I3P1

						   !if (I1>I2) cycle
						   !if (I1==I2 .and. I23>I13) cycle

						   INX =INX +1
						   PowerTable(INX,1) = I1
						   PowerTable(INX,2) = I2
						   PowerTable(INX,3) = I12
						   PowerTable(INX,4) = I3
						   PowerTable(INX,5) = I13
						   PowerTable(INX,6) = I23

						   !write (*,*) I1, I2, I12, I3, I13, I23

	70                  CONTINUE
	69                CONTINUE
	68              CONTINUE
	67            CONTINUE
	66          CONTINUE
	65        CONTINUE
		NumTerm = INX
	endif

! REMOVE THIS !
	if (Ordering == 2) then  ! Use my ordering
		NumTerm = 0
		do om = 0, Omega, 1
			do qi = 0, Omega, 1
				do pi = 0, Omega, 1
					do ni = 0, Omega, 1
						do mi = 0, Omega, 1
							do li = 0, Omega, 1
								do ki = 0, Omega, 1
									if (ki + li + mi + ni + pi + qi == om) then
										!if (ki + li + mi + ni + pi + qi <= Omega) then
										!if (ki>li) cycle
										!if (ki==li .and. qi>pi) cycle

	!									if (om > 3 .and. qi > 0) cycle  ! Restrict r23 to max of 1
	!									if (om > 3 .and. pi > 0) cycle  ! Restrict r13 to max of 3
	!									if (om > 2 .and. ni > 0) cycle  ! Restrict r3 to max of 0
	!									if (ki >= li .and. qi > pi) cycle  ! Check if r1>r2 and r23 > r13 (powers)

										NumTerm = NumTerm + 1
										PowerTable(NumTerm,1) = ki
										PowerTable(NumTerm,2) = li
										PowerTable(NumTerm,3) = ni
										PowerTable(NumTerm,4) = mi
										PowerTable(NumTerm,5) = qi
										PowerTable(NumTerm,6) = pi

										write (*,*) ki, li, ni, mi, qi, pi!, ki + li + mi + ni + pi + qi
									endif
								enddo
							enddo
						enddo
					enddo
				enddo
			enddo
		enddo
	endif

	do i = 1, NumTerm, 1
		! Increase r1 by l1 and r2 by l2/l3
		PowerTable(i,1) = PowerTable(i,1) + l1
		PowerTable(i,2) = PowerTable(i,2) + l2
		write (*,'(i3,a,6i2)') i, ":", PowerTable(i,1), PowerTable(i,2), PowerTable(i,3), PowerTable(i,4), PowerTable(i,5), PowerTable(i,6)
	end do

	return
end

subroutine GenCoeffTable(CoeffTable, PowerTablei, i, PowerTablej, j, NumTerms, Alphai, Alphaj, Betai, Betaj, &
						 Gammai, Gammaj, l1r, l2r, l3r)
	implicit none
	real*16, dimension(34) :: CoeffTable
	integer i, j, NumTerms
	integer, dimension(NumTerms,6) :: PowerTablei
	integer, dimension(NumTerms,6) :: PowerTablej
	real*16 Alphai, Alphaj, Betai, Betaj, Gammai, Gammaj
	integer kj, lj, mj, nj, pj, qj, l1r, l2r, l3r

	! We do not really have to do this, but it makes reading the assignments later
	!  easier to follow instead of assigning index numbers to each of k, l, m, etc.
	! In the power tables, we have the following indices:
	!  k = 1, l = 2, m = 3, n = 4, p = 5, q = 6
	kj = PowerTablej(j,1)
	lj = PowerTablej(j,2)
	mj = PowerTablej(j,3)
	nj = PowerTablej(j,4)
	pj = PowerTablej(j,5)
	qj = PowerTablej(j,6)

	CoeffTable(1) = -Alphaj*Alphaj - Betaj*Betaj - Gammaj*Gammaj
	CoeffTable(2) = -kj - kj*kj - kj*mj - kj*pj + l1r + l1r*l1r
	CoeffTable(3) = -2.0q0*mj - kj*mj - 2.0q0*mj*mj - mj*lj - mj*qj - mj*pj
	CoeffTable(4) = -2.0q0
	CoeffTable(5) = -lj - mj*lj - lj*lj - lj*qj + l2r + l2r*l2r
	CoeffTable(6) = mj*lj
	CoeffTable(7) = kj*mj
	CoeffTable(8) = -2.0q0*qj - mj*qj - lj*qj - 2.0q0*qj*qj - qj*nj - qj*pj
	CoeffTable(9) = 2.0q0
	CoeffTable(10) = -nj - qj*nj - nj*nj - nj*pj + l3r + l3r*l3r
	CoeffTable(11) = qj*nj
	CoeffTable(12) = lj*qj
	CoeffTable(13) = -2.0q0*pj - kj*pj - mj*pj - qj*pj - nj*pj - 2.0q0*pj*pj
	CoeffTable(14) = qj*pj
	CoeffTable(15) = mj*pj
	CoeffTable(16) = nj*pj
	CoeffTable(17) = kj*pj
	CoeffTable(18) = -2.0q0
	CoeffTable(19) = mj*qj
	CoeffTable(20) = mj*Alphaj
	CoeffTable(21) = -mj*Alphaj
	CoeffTable(22) = pj*Alphaj
	CoeffTable(23) = -pj*Alphaj
	CoeffTable(24) = 2.0q0 + 2.0q0*Alphaj + 2.0q0*kj*Alphaj + mj*Alphaj + pj*Alphaj
	CoeffTable(25) = -mj*Betaj
	CoeffTable(26) = mj*Betaj
	CoeffTable(27) = qj*Betaj
	CoeffTable(28) = -qj*Betaj
	CoeffTable(29) = -2.0q0 + 2.0q0*Betaj + mj*Betaj + 2.0q0*lj*Betaj + qj*Betaj
	CoeffTable(30) = -qj*Gammaj
	CoeffTable(31) = qj*Gammaj
	CoeffTable(32) = -pj*Gammaj
	CoeffTable(33) = pj*Gammaj
	CoeffTable(34) = -2.0q0 + 2.0q0*Gammaj + qj*Gammaj + 2.0q0*nj*Gammaj + pj*Gammaj

	return
end

! Note that the naming is a little different here.  I used rPowerTable here as the parameter name but
!  use rPowers elsewhere.
subroutine GenRPowerTable(rPowerTable, PowerTablei, i, PowerTablej, j, NumTerms)
	integer, dimension(34,6) :: rPowerTable
	integer, dimension(NumTerms,6) :: PowerTablei
	integer, dimension(NumTerms,6) :: PowerTablej
	integer i, j, NumTerms
	integer ki, li, mi, ni, pi, qi, kj, lj, mj, nj, pj, qj
	integer, dimension(34,6) :: rPowers
	integer n

	! This is a 2-dimensional array expressed in a contiguous 1-D array.  I see no way to assign a 2-D
	!  array in Fortran (at least in G95).  I found several examples that would not compile with G95.
	!  These are taken directly from the table in rpowers.pdf.  Each power is covered by 2 lines, since
	!  there are 34 terms to take care of.  These can be read as going down the columns of rpowers.pdf.
	data rPowers / 0,-2, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 1,-1, 1,-1,-1, 2, 0, 0, 0, 0, 0, 0, 2, 0, 0, &  ! r1
				   0, 0, 0, 0,-2,-2, 2, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-1, 1, 1,-1,-1, 2, 0, 0, 0, 0, &  ! r2
				   0, 0,-2,-1, 0,-2,-2, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0,-2,-2,-2, 0, 0, 0,-2,-2, 0, 0, 0, 0, 0, 0, 0, 0, &  ! r12
				   0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2, 2, 0, 0, 0,-2, 2, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 2, 0,-1, 1,-1, 1,-1, &  ! r3
				   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2,-2,-2,-2,-1, 2, 0, 0,-2,-2, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2, 0, &  ! r31
				   0, 0, 0, 0, 0, 0, 0,-2,-1, 0,-2,-2, 0,-2, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0,-2,-2, 0,-2,-2, 0, 0, 0 /   ! r23
				  !1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34

	! We do not really have to do this, but it makes reading the assignments later
	!  easier to follow instead of assigning index numbers to each of k, l, m, etc.
	! In the power tables, we have the following indices:
	!  k = 1, l = 2, m = 3, n = 4, p = 5, q = 6
	ki = PowerTablei(i,1)
	li = PowerTablei(i,2)
	mi = PowerTablei(i,3)
	ni = PowerTablei(i,4)
	pi = PowerTablei(i,5)
	qi = PowerTablei(i,6)
	kj = PowerTablej(j,1)
	lj = PowerTablej(j,2)
	mj = PowerTablej(j,3)
	nj = PowerTablej(j,4)
	pj = PowerTablej(j,5)
	qj = PowerTablej(j,6)

	! Combine powers of like r's from f_i, f_j and the powers in rpowers.pdf.
	do n = 1, 34, 1
		rPowerTable(n,1) = ki + kj + rPowers(n,1)
		rPowerTable(n,2) = li + lj + rPowers(n,2)
		rPowerTable(n,3) = mi + mj + rPowers(n,3)
		rPowerTable(n,4) = ni + nj + rPowers(n,4)
		rPowerTable(n,5) = pi + pj + rPowers(n,5)
		rPowerTable(n,6) = qi + qj + rPowers(n,6)
	enddo

	return
end

subroutine GenCoeffTableSH(CoeffTable, PowerTablei, i, PowerTablej, j, NumTerms, Alphai, Alphaj, Betai, Betaj, Gammai, Gammaj, l1r, l2r, l3r)
	implicit none
	real*16, dimension(6) :: CoeffTable
	integer i, j, NumTerms
	integer, dimension(NumTerms,6) :: PowerTablei
	integer, dimension(NumTerms,6) :: PowerTablej
	real*16 Alphai, Alphaj, Betai, Betaj, Gammai, Gammaj
	integer kj, lj, mj, nj, pj, qj, l1r, l2r, l3r

	kj = PowerTablej(j,1)
	lj = PowerTablej(j,2)
	mj = PowerTablej(j,3)
	nj = PowerTablej(j,4)
	pj = PowerTablej(j,5)
	qj = PowerTablej(j,6)

	CoeffTable(1) = mj
	CoeffTable(2) = pj
	CoeffTable(3) = qj
	CoeffTable(4) = mj
	CoeffTable(5) = pj
	CoeffTable(6) = qj

	return
end

! Note that the naming is a little different here.  I used rPowerTable here as the parameter name but
!  use rPowers elsewhere.
subroutine GenRPowerTableSH(rPowerTable, PowerTablei, i, PowerTablej, j, NumTerms)
	integer, dimension(6,6) :: rPowerTable
	integer, dimension(NumTerms,6) :: PowerTablei
	integer, dimension(NumTerms,6) :: PowerTablej
	integer i, j, NumTerms
	integer ki, li, mi, ni, pi, qi, kj, lj, mj, nj, pj, qj
	integer, dimension(6,6) :: rPowers
	integer n

	data rPowers / -1,-1, 0, 1, 1, 0, &  ! r1
				    1, 0,-1,-1, 0, 1, &  ! r2
				   -2, 0, 0,-2, 0, 0, &  ! r12
				    0, 1, 1, 0,-1,-1, &  ! r3
				    0,-2, 0, 0,-2, 0, &  ! r31
				    0, 0,-2, 0, 0,-2 /   ! r23
				   !1, 2, 3, 4, 5, 6

	! In the power tables, we have the following indices:
	!  k = 1, l = 2, m = 3, n = 4, p = 5, q = 6
	ki = PowerTablei(i,1)
	li = PowerTablei(i,2)
	mi = PowerTablei(i,3)
	ni = PowerTablei(i,4)
	pi = PowerTablei(i,5)
	qi = PowerTablei(i,6)
	kj = PowerTablej(j,1)
	lj = PowerTablej(j,2)
	mj = PowerTablej(j,3)
	nj = PowerTablej(j,4)
	pj = PowerTablej(j,5)
	qj = PowerTablej(j,6)

	! Combine powers of like r's from f_i, f_j and the powers in rpowers.pdf.
	do n = 1, 6, 1
		rPowerTable(n,1) = ki + kj + rPowers(n,1)
		rPowerTable(n,2) = li + lj + rPowers(n,2)
		rPowerTable(n,3) = mi + mj + rPowers(n,3)
		rPowerTable(n,4) = ni + nj + rPowers(n,4)
		rPowerTable(n,5) = pi + pj + rPowers(n,5)
		rPowerTable(n,6) = qi + qj + rPowers(n,6)
	enddo

	return
end

! Sets all WMatrix elements equal to 0.
! TODO: Is there an equivalent to the C++ memset?
subroutine ClearWMatrix(WMatrix)
	use WLimits
	implicit none
	real*16, dimension(lmin:lmax, mmin:mmax, nmin:nmax, 6) :: WMatrix
	integer l, m, n

	do l = lmin, lmax, 1
		do m = mmin, mmax, 1
			do n = nmin, nmax, 1
				WMatrix(l, m, n, 1) = 0.0q0
				WMatrix(l, m, n, 2) = 0.0q0
				WMatrix(l, m, n, 3) = 0.0q0
				WMatrix(l, m, n, 4) = 0.0q0
				WMatrix(l, m, n, 5) = 0.0q0
				WMatrix(l, m, n, 6) = 0.0q0
			end do
		end do
	end do

	return
end

! This precalculates the W matrices, since the W function is computationally expensive.
subroutine CalcWMatrices(Omega, WMatrix, Alpha, Beta, Gamma, pmax)
	use WLimits
	implicit none
	integer Omega
	real*16, dimension(lmin:lmax, mmin:mmax, nmin:nmax, 6) :: WMatrix
	real*16 Alpha, Beta, Gamma
	integer pmax
	real*16 W
	integer l, m, n, k

	!$omp parallel do shared(lmin,lmax,mmin,mmax,nmin,nmax,pmax,Alpha,Beta,Gamma,WMatrix) private(m,n) schedule(dynamic,5)
	do l = lmin, lmax, 1
		write (*,*) "wmatrix l:", l, "/", lmax
		!write (*,"(a)",advance='no') '.'
		do m = mmin, mmax, 1
			do n = nmin, nmax, 1
				!if (l + m + n + 2 >= 0 .and. l + m + 1 >= 0) then  ! TODO: Needed?
					do k = 1, 6, 1
						! Note that the ordering of alpha, beta and gamma is the same as the ordering
						!  in the terms of equation (6) of Drake/Yan '95.
						if (WMatrix(l, m, n, k) /= 0.0q0) then
							select case(k)
								case (1)
									WMatrix(l, m, n, 1) = W(l, m, n, Alpha, Beta, Gamma, pmax)
								case (2)
									WMatrix(l, m, n, 2) = W(l, m, n, Alpha, Gamma, Beta, pmax)
								case (3)
									WMatrix(l, m, n, 3) = W(l, m, n, Beta, Alpha, Gamma, pmax)
								case (4)
									WMatrix(l, m, n, 4) = W(l, m, n, Beta, Gamma, Alpha, pmax)
								case (5)
									WMatrix(l, m, n, 5) = W(l, m, n, Gamma, Alpha, Beta, pmax)
								case (6)
									WMatrix(l, m, n, 6) = W(l, m, n, Gamma, Beta, Alpha, pmax)
							end select  ! There is no need for a default case.
						end if
					end do
				!end if
			end do
		end do
	end do
	write (*,*)

	return
end

subroutine CalcMatricesSub(RunCalc, UsePreCalc, CalcSH, PowerTablei, PowerTablej, WMatrix, PhiPhi, Phi2HPhi, NumTerms, OffI, OffJ, Alphai, Alphaj, Betai, Betaj, &
							Gammai, Gammaj, l1l, l2l, l3l, m1l, m2l, m3l, l1r, l2r, l3r, m1r, m2r, m3r, M12max, M23max, M31max, pmax, IsTriplet, Method)
	use WLimits
	implicit none
	real*16, dimension(34) :: CoeffTable
	real*16, dimension(6) :: CoeffTableSH
	integer, dimension(34,6) :: rPowers
	integer, dimension(6,6) :: rPowersSH
	integer, dimension(6) :: sldata, srdata
	integer, dimension(NumTerms,6) :: PowerTablei, PowerTablej
	real*16, dimension(lmin:lmax, mmin:mmax, nmin:nmax, 6) :: WMatrix
	real*16, dimension(NumTerms*2,NumTerms*2) :: PhiPhi, Phi2HPhi
	integer M12max, M23max, M31max, pmax, NumTerms, OffI, OffJ, Method
	logical RunCalc, UsePreCalc, CalcSH
	integer l1l, l2l, l3l, m1l, m2l, m3l, l1r, l2r, l3r, m1r, m2r, m3r
	real*16 Alphai, Alphaj, Betai, Betaj, Gammai, Gammaj
	real*16 HylleraasIntegral, HylleraasIntegralGeneral, CP12Angular
	real*16 Sum, SumPiHalf
	real*16 PI
	integer i, j, n, IsTriplet
	real*16 RemoveMe, SHPart, SHPartPiHalf
	real*16 PiHalf
	data sldata / 2, 3, 3, 1, 1, 2 /
	data srdata / 1, 1, 2, 2, 3, 3 /

	! 1/(2pi) calculated in Mathematica
	PiHalf = 0.1591549430918953357688837633725143620345_16

	!call omp_set_num_threads(1)

	!$omp parallel do shared(NumTerms) private(j,Sum,SHPart,RemoveMe,CoeffTable,CoeffTableSH,rPowers,rPowersSH) schedule(dynamic,10)
	do i = 1, NumTerms, 1
		write (*,*) i
		!write (*,"(i7)",advance='no') i
		!do j = 1, i, 1
		do j = 1, NumTerms, 1
			Sum = 0.0q0
			call GenCoeffTable(CoeffTable, PowerTablei, i, PowerTablej, j, NumTerms, &
								Alphai, Alphaj, Betai, Betaj, Gammai, Gammaj, l1r, l2r, l3r)
			call GenRPowerTable(rPowers, PowerTablei, i, PowerTablej, j, NumTerms)
			call GenCoeffTableSH(CoeffTableSH, PowerTablei, i, PowerTablej, j, NumTerms, &
								Alphai, Alphaj, Betai, Betaj, Gammai, Gammaj, l1r, l2r, l3r)
			call GenRPowerTableSH(rPowersSH, PowerTablei, i, PowerTablej, j, NumTerms)

			do n = 1, 34, 1
            !do n = 4, 4, 1
				if (CoeffTable(n) /= 0.0q0) then
					! Note that the order of the parameters to the HylleraasIntegral function looks strange. This is
					!  because the Drake and Yan paper (1995) expressed the integrand in the following order:
					!  r1 r2 r3 r12 r23 r31, and the Armour and Humberston article has them in this order:
					!  r1 r2 r12 r3 r13 r23 (equation 3.14). So we have to swap the third and fourth parameters, along
					!  with the fifth and sixth.

					! Do the PhiPhi inner product
					RemoveMe = HylleraasIntegralGeneral(RunCalc, UsePreCalc, rPowers(n,1), rPowers(n,2), rPowers(n,4), rPowers(n,3), rPowers(n,6), &
												  rPowers(n,5), Alphai+Alphaj, Betai+Betaj, Gammai+Gammaj, &
												  l1l, l2l, l3l, m1l, m2l, m3l, l1r, l2r, l3r, m1r, m2r, m3r, M12max, M23max, M31max, pmax, 0, 0, 0, WMatrix, Method)
					Sum = Sum + RemoveMe * CoeffTable(n)
				end if
			enddo

			! TODO: Check for whether these terms should be 0 before calculating.
			! Now we need to take care of the terms that act on the spherical harmonics. There are 6 terms total from this.
			if (CalcSH == .true.) then
				SHPart = 0.0q0
				do n = 1, 6, 1
					if (CoeffTableSH(n) /= 0.0q0) then
						RemoveMe = HylleraasIntegralGeneral(RunCalc, UsePreCalc, rPowersSH(n,1), rPowersSH(n,2), rPowersSH(n,4), rPowersSH(n,3), rPowersSH(n,6), &
															  rPowersSH(n,5), Alphai+Alphaj, Betai+Betaj, Gammai+Gammaj, &
															  l1l, l2l, l3l, m1l, m2l, m3l, l1r, l2r, l3r, m1r, m2r, m3r, M12max, M23max, M31max, pmax, n, sldata(n), srdata(n), WMatrix, Method)
						SHPart = SHPart + RemoveMe * CoeffTableSH(n)
					end if
				end do
				!SumPiHalf = Sum * PiHalf  ! TODO: Remove this.
				!SHPartPiHalf = SHPart * PiHalf  ! TODO: Remove this.
				Sum = Sum + SHPart * 2  ! The other terms have a 1/2 in front, and we are calculating 2H
			end if

			if (IsTriplet == 0) then
				Phi2HPhi(i+OffI,j+OffJ) = Phi2HPhi(i+OffI,j+OffJ) + Sum * PiHalf
			else
				Phi2HPhi(i+OffI,j+OffJ) = Phi2HPhi(i+OffI,j+OffJ) - Sum * PiHalf
			end if

			! Do the PhiPhi inner product
			RemoveMe = HylleraasIntegralGeneral(RunCalc, UsePreCalc, PowerTablei(i,1)+PowerTablej(j,1), PowerTablei(i,2)+PowerTablej(j,2), &
						PowerTablei(i,4)+PowerTablej(j,4), PowerTablei(i,3)+PowerTablej(j,3), &
						PowerTablei(i,6)+PowerTablej(j,6), PowerTablei(i,5)+PowerTablej(j,5), Alphai+Alphaj, Betai+Betaj, Gammai+Gammaj, &
						l1l, l2l, l3l, m1l, m2l, m3l, l1r, l2r, l3r, m1r, m2r, m3r, M12max, M23max, M31max, pmax, 0, 0, 0, WMatrix, Method)
			if (IsTriplet == 0) then
				PhiPhi(i+OffI,j+OffJ) = PhiPhi(i+OffI,j+OffJ) + RemoveMe * PiHalf
			else
				PhiPhi(i+OffI,j+OffJ) = PhiPhi(i+OffI,j+OffJ) - RemoveMe * PiHalf
			endif
		enddo
	enddo

	write (*,*)
	return
end

subroutine CalcMatrices(RunCalc, UsePreCalc, PowerTabler1i, PowerTabler1j, PowerTabler2i, PowerTabler2j, WMatrix, PhiPhi, Phi2HPhi, NumTerms, &
						Alphai, Alphaj, Betai, Betaj, Gammai, Gammaj, LValue, M12max, M23max, M31max, pmax, IsTriplet, Method, Exchanged)
	use WLimits
	implicit none
	integer, dimension(NumTerms/2,6) :: PowerTabler1i, PowerTabler1j, PowerTabler2i, PowerTabler2j
	real*16, dimension(lmin:lmax, mmin:mmax, nmin:nmax, 6) :: WMatrix
	real*16, dimension(NumTerms,NumTerms) :: PhiPhi, Phi2HPhi
	integer NumTerms, LValue, M12max, M23max, M31max, pmax, IsTriplet, Method
	logical RunCalc, UsePreCalc
	real*16 Alphai, Alphaj, Betai, Betaj, Gammai, Gammaj
	logical Exchanged

	! 1-1 matrix elements
	write (*,*) "Calculating 1-1 matrix elements"
	call CalcMatricesSub(RunCalc, UsePreCalc, .false., PowerTabler1i, PowerTabler1j, WMatrix, PhiPhi, Phi2HPhi, NumTerms/2, 0, 0, Alphai, Alphaj, Betai, Betaj, &
						Gammai, Gammaj, LValue, 0, 0, 0, 0, 0, LValue, 0, 0, 0, 0, 0, M12max, M23max, M31max, pmax, IsTriplet, Method)

    if (LValue == 0) then  ! S-wave only has one symmetry
        return
    end if

	! 2-2 matrix elements
	write (*,*) "Calculating 2-2 matrix elements"
	if (Exchanged == .true.) then  ! Swap l2r and l3r
		call CalcMatricesSub(RunCalc, UsePreCalc, .true., PowerTabler2i, PowerTabler2j, WMatrix, PhiPhi, Phi2HPhi, NumTerms/2, NumTerms/2, NumTerms/2, Alphai, Alphaj, Betai, Betaj, &
							Gammai, Gammaj, 0, LValue, 0, 0, 0, 0, 0, 0, LValue, 0, 0, 0, M12max, M23max, M31max, pmax, IsTriplet, Method)
	else
		call CalcMatricesSub(RunCalc, UsePreCalc, .false., PowerTabler2i, PowerTabler2j, WMatrix, PhiPhi, Phi2HPhi, NumTerms/2, NumTerms/2, NumTerms/2, Alphai, Alphaj, Betai, Betaj, &
							Gammai, Gammaj, 0, LValue, 0, 0, 0, 0, 0, LValue, 0, 0, 0, 0, M12max, M23max, M31max, pmax, IsTriplet, Method)
	end if

	! 1-2 matrix elements
	write (*,*) "Calculating 1-2 matrix elements"
	if (Exchanged == .true.) then  ! Swap l2r and l3r
		call CalcMatricesSub(RunCalc, UsePreCalc, .true., PowerTabler1i, PowerTabler2j, WMatrix, PhiPhi, Phi2HPhi, NumTerms/2, 0, NumTerms/2, Alphai, Alphaj, Betai, Betaj, &
							Gammai, Gammaj, LValue, 0, 0, 0, 0, 0, 0, 0, LValue, 0, 0, 0, M12max, M23max, M31max, pmax, IsTriplet, Method)
	else
		call CalcMatricesSub(RunCalc, UsePreCalc, .true.,  PowerTabler1i, PowerTabler2j, WMatrix, PhiPhi, Phi2HPhi, NumTerms/2, 0, NumTerms/2, Alphai, Alphaj, Betai, Betaj, &
							Gammai, Gammaj, LValue, 0, 0, 0, 0, 0, 0, LValue, 0, 0, 0, 0, M12max, M23max, M31max, pmax, IsTriplet, Method)
	end if

	! 2-1 matrix elements
	write (*,*) "Calculating 2-1 matrix elements"
	call CalcMatricesSub(RunCalc, UsePreCalc, .true., PowerTabler2i, PowerTabler1j, WMatrix, PhiPhi, Phi2HPhi, NumTerms/2, NumTerms/2, 0, Alphai, Alphaj, Betai, Betaj, &
						Gammai, Gammaj, 0, LValue, 0, 0, 0, 0, LValue, 0, 0, 0, 0, 0, M12max, M23max, M31max, pmax, IsTriplet, Method)

	return
end

!real*16 function ederivwrapper(Omega, NumTerms, M12max, pmax, IsTriplet, Ordering, Method, EigenRoutine, Alpha, Beta, Gamma, Iter, Energy)
!	implicit none
!	real*16, allocatable, dimension(:,:) :: PhiPhi, Phi2HPhi
!	real*16, allocatable, dimension(:,:) :: PhiPhi8, Phi2HPhi8
!	real*16, allocatable, dimension(:) :: Energies, Workspace
!	integer Omega, NumTerms, M12max, pmax, IsTriplet, Ordering, Method, EigenRoutine, Info, Iter
!	real*16 Alpha, Beta, Gamma, Energy
!	real*16 EDeriv
!
!	allocate(PhiPhi(NumTerms,NumTerms))
!	allocate(Phi2HPhi(NumTerms,NumTerms))
!	allocate(Energies(NumTerms))
!	allocate(Workspace(3*NumTerms-1))  ! Not really needed for Newdsygv?
!
!	Energies = 0.0q0
!
!	Energy = EDeriv(Omega, NumTerms, Alpha, Beta, Gamma, 1, PhiPhi, Phi2HPhi, PhiPhi8, Phi2HPhi8, &
!					M12max, pmax, 0, IsTriplet, Ordering, Method, EigenRoutine)
!
!	Phi2HPhi = Phi2HPhi / 2.0q0  ! TODO: Remove?
!	if (EigenRoutine == 1) then
!		allocate(PhiPhi8(NumTerms,NumTerms))
!		allocate(Phi2HPhi8(NumTerms,NumTerms))
!		PhiPhi8 = PhiPhi
!		Phi2HPhi8 = Phi2HPhi
!		call dsygv(1, 'N', 'L', NumTerms, Phi2HPhi8, NumTerms, PhiPhi8, NumTerms, Energies, Workspace, 3*NumTerms-1, Info)
!		! TODO: Test Info!
!		write (*,*) "Info: ", Info
!		deallocate(PhiPhi8, Phi2HPhi8)
!	else  ! EigenRoutine == 2
!		call Newdsygv(1, 'N', 'L', NumTerms, Phi2HPhi, NumTerms, PhiPhi, NumTerms, Energies, Workspace, 3*NumTerms-1, Info)
!	end if
!
!	Energy = Energies(1)
!
!	write (*,"(A,3F8.5)") "Alpha, Beta, Gamma: ", Alpha, Beta, Gamma
!	write (10,"(I4,A,3F8.5)") Iter, "; Alpha, Beta, Gamma: ", Alpha, Beta, Gamma  ! I dislike hardcoding the 10, but it's easier at this point.
!	write (10,"(A,F16.12)") "Energy: ", Energy
!	write (10,*)
!	flush(10)
!
!	deallocate(PhiPhi, Phi2HPhi, Energies, Workspace)
!	ederivwrapper = 1.0q0
!end

! Calculates the derivative of the energy function.  If the Deriv parameter = 0, then it exits early with
!  PhiPhi and Phi2HPhi filled.
real*16 function EDeriv(Omega, NumTerms, Alpha, Beta, Gamma, EigenNum, PhiPhi, Phi2HPhi, &
						M12max, M23max, M31max, pmax, Deriv, IsTriplet, Ordering, Method, EigenRoutine, LValue)
	use WLimits
	!use nag_nsym_gen_eig
	implicit none
	real*16 Alpha, Beta, Gamma, Alphai, Alphaj, Betai, Betaj, Gammai, Gammaj, e
	integer, allocatable, dimension(:,:) :: PowerTabler1i, PowerTabler1j, PowerTabler2i, PowerTabler2j  ! This contains the powers of the Hylleraas terms.
	integer CalcPowerTableSize, PreCalcGammaSize  ! Function definitions
	integer EigenNum  ! This is the energy eigenvalue number we are dealing with.
	integer Omega, NumTerms, M12max, M23max, M31max, pmax, Deriv, IsTriplet, Method, Ordering, i, j, n, TempJ, Info, EigenRoutine, LValue
	real*16, dimension(NumTerms,NumTerms) :: PhiPhi, Phi2HPhi
	real*16, allocatable, dimension(:,:,:,:) :: WMatrix
	real*16 PhiHPhiSum, PhiPhiSum, Norm
	integer l, m, k
	logical UsePreCalc

	UsePreCalc = .true.
	!UsePreCalc = .false.

	allocate(PowerTabler1i(NumTerms/2, 6))
	call GenOmegaPowerTable(Omega, PowerTabler1i, NumTerms/2, Ordering, LValue, 0)
	allocate(PowerTabler1j(NumTerms/2, 6))
	call GenOmegaPowerTable(Omega, PowerTabler1j, NumTerms/2, Ordering, LValue, 0)
	allocate(PowerTabler2i(NumTerms/2, 6))
	call GenOmegaPowerTable(Omega, PowerTabler2i, NumTerms/2, Ordering, 0, LValue)
	allocate(PowerTabler2j(NumTerms/2, 6))
	call GenOmegaPowerTable(Omega, PowerTabler2j, NumTerms/2, Ordering, 0, LValue)

	! Set this matrix equal to 0 because of the Phi2HPhi(i,j) = Phi2HPhi(i,j) + Sum line in CalcMatrices
	!  (and equivalent for PhiPhi).  This allows us to use the same subroutine for the direct and exchange elements.
	Phi2HPhi = 0.0_16
	PhiPhi = 0.0_16

	! For the direct integration, f_i and f_j have beta and gamma in the same place.
	Alphai = Alpha; Alphaj = Alpha; Betai = Beta; Betaj = Beta; Gammai = Gamma; Gammaj = Gamma;

	! Precalculate the W matrices.  lmin, lmax, etc. are the upper and lower limits for each of the three
	!  parameters to the W function (in HylleraasIntegral.f90).  These limits were determined by me, and
	!  there is a document detailing these limits (W Function Limits.pdf).
	!!lmin = 1;  lmax = 2*(Omega+LValue) + 4*max(linf,minf) + 4;
	!!mmin = -1 - 3*max(linf,minf);  mmax = 2*(Omega+LValue) + 3*max(linf,minf) + 3;
	!!nmin = -2 - 4*max(linf,minf);  nmax = 2*(Omega+LValue) + 2;
	lmin = 1;  lmax = min(M12max, min(M23max, M31max))*2 + LValue*2 + Omega*2 + 5;
	mmin = -1;  mmax = 2*(LValue+Omega) + 3;
	nmin = -(min(M12max, min(M23max, M31max))*3 + 2*Omega);  nmax = 2*(Omega+LValue) + 2;
	allocate(WMatrix(lmin:lmax, mmin:mmax, nmin:nmax, 6))

    call GenFactorial(400)

	write (*,*)
	write (*,*) "Starting calculation of direct-direct terms"
	if (UsePreCalc == .true.) then
		write (*,*) "Precomputing direct-direct W matrix"
		call ClearWMatrix(WMatrix)
		! First call it with RunCalc = false to determine what W functions to calculate.
		call CalcMatrices(.false., .true., PowerTabler1i, PowerTabler1j, PowerTabler2i, PowerTabler2j, WMatrix, PhiPhi, Phi2HPhi, NumTerms, &
							Alphai, Alphaj, Betai, Betaj, Gammai, Gammaj, LValue, M12max, M23max, M31max, pmax, IsTriplet, Method, .false.)
		! Then precalculate all W functions that are used.
		call CalcWMatrices(Omega, WMatrix, Alphai+Alphaj, Betai+Betaj, Gammai+Gammaj, pmax)
		write (*,*) "Finished precomputing direct-direct W matrix"
	end if
	! This subroutine calculates the PhiPhi and Phi2HPhi matrices for the direct integrals.
	call CalcMatrices(.true., UsePreCalc, PowerTabler1i, PowerTabler1j, PowerTabler2i, PowerTabler2j, WMatrix, PhiPhi, Phi2HPhi, NumTerms, &
						Alphai, Alphaj, Betai, Betaj, Gammai, Gammaj, LValue, M12max, M23max, M31max, pmax, IsTriplet, Method, .false.)
	if (IsTriplet == 1) then  ! Only do this for the triplet case.
		PhiPhi = -PhiPhi
		Phi2HPhi = -Phi2HPhi
    endif

	!goto 300  ! Skip permutation for testing

	! Now we have to do the same thing for the exchange operator, P_23.  We operate it on f_j by switching the
	!  2<->3 powers in PowerTablej and recalculating, along with beta and gamma.
	do n = 1, NumTerms/2, 1
		! This exchanges r2 and r3.
		TempJ = PowerTabler1j(n,2)
		PowerTabler1j(n,2) = PowerTabler1j(n,4)  ! l with n
		PowerTabler1j(n,4) = TempJ
		TempJ = PowerTabler2j(n,2)
		PowerTabler2j(n,2) = PowerTabler2j(n,4)  ! l with n
		PowerTabler2j(n,4) = TempJ
		! This exchanges r12 and r13 (r31).
		TempJ = PowerTabler1j(n,3)
		PowerTabler1j(n,3) = PowerTabler1j(n,5)  ! m with p
		PowerTabler1j(n,5) = TempJ
		TempJ = PowerTabler2j(n,3)
		PowerTabler2j(n,3) = PowerTabler2j(n,5)  ! m with p
		PowerTabler2j(n,5) = TempJ
	enddo
	! Swap beta and gamma for f_j.
	Alphai = Alpha; Alphaj = Alpha; Betaj = Gamma; Betai = Beta; Gammaj = Beta; Gammai = Gamma;

	!i = 0
	!do l = lmin, lmax, 1
	!	do m = mmin, mmax, 1
	!		do n = nmin, nmax, 1
	!			do k = 1, 6, 1
	!				if (WMatrix(l, m, n, k) == 1000.0q0) i = i + 1
	!			end do
	!		end do
	!	end do
	!end do
	!
	!write (*,*) "Total WMatrix used: ", i
	!write (*,*) "Total WMatrix size: ", (lmax-lmin+1) * (mmax-mmin+1) * (nmax-nmin+1) * 6
	!stop

	! We have to recalculate the W matrices again, since we had the swaps.
	write (*,*) "Starting calculation of direct-exchange terms"
	if (UsePreCalc == .true.) then
		write (*,*) "Precomputing direct-exchange W matrix"
		call ClearWMatrix(WMatrix)
		! The W functions to calculate are not, in general, the same as that for the direct-direct case.
		call CalcMatrices(.false., UsePreCalc, PowerTabler1i, PowerTabler1j, PowerTabler2i, PowerTabler2j, WMatrix, PhiPhi, Phi2HPhi, NumTerms, &
						Alphai, Alphaj, Betai, Betaj, Gammai, Gammaj, LValue, M12max, M23max, M31max, pmax, IsTriplet, Method, .true.)
		call CalcWMatrices(Omega, WMatrix, Alphai+Alphaj, Betai+Betaj, Gammai+Gammaj, pmax)
		write (*,*) "Finished precomputing direct-exchange W matrix"
	end if
	! Calculate the matrix elements for the exchanged terms, and add them to the direct terms.
	call CalcMatrices(.true., UsePreCalc, PowerTabler1i, PowerTabler1j, PowerTabler2i, PowerTabler2j, WMatrix, PhiPhi, Phi2HPhi, NumTerms, &
						Alphai, Alphaj, Betai, Betaj, Gammai, Gammaj, LValue, M12max, M23max, M31max, pmax, IsTriplet, Method, .true.)

	! If we are not optimizing the nonlinear parameters, return since matrices have been calculated.
300	EDeriv = 1.0q0  ! Random value, return not used in this case.

	! Clean up memory before exiting
	deallocate(WMatrix)
	deallocate(PowerTabler1i)
	deallocate(PowerTabler1j)
	deallocate(PowerTabler2i)
	deallocate(PowerTabler2j)

	return
end

max_stars_repo_name: DentonW/Ps-H-Scattering

max_issues_repo_name: DentonW/Ps-H-Scattering

max_forks_repo_name: DentonW/Ps-H-Scattering


text:
C++***************************************************************
C Program spec_extract
C To extract a spectrum from a 2D image
C Linear interpolation for slice extraction
C Extraction of a spectrum:
C  - from a 2D image containing only one spectrum if rotation is needed
C  - or containing more spectra if horizontal case(please run "resample1" before)')
C
C JLP
C Version of 17-12-98
C--**************************************************************
	PROGRAM SPEC_EXTRACT
	CHARACTER REP*1,IN_NAME*40,IN_COMMENTS*80
	CHARACTER OUT_NAME*40,OUT_COMMENTS*80
	INTEGER*4 MADRID(1),PNTR_1,PNTR_2
	INTEGER*4 IX1,IX2,IY1,IY2,NX_2,NY_2
	INTEGER*4 NX_1,NY_1
	COMMON /VMR/MADRID

10	FORMAT(A)

	CALL JLP_BEGIN

	WRITE(6,22)
22      FORMAT(' Program spec_extract  Version of 17-12-98',/,
     1         ' Extraction from a 2D image containing only one spectrum',
     1         ' if rotation is applied',/,
     1         ' or more spectra if horizontal (please run "resample1" before)')
        WRITE(6,*) ' JLP99: Inversion of wavelengths for speckle data'

C Inquires about the format of the files :
	CALL JLP_INQUIFMT

C************************************************************
C Loading input image file
C************************************************************
        WRITE(6,*) 'Input file: '
        READ(5,10) IN_NAME
	CALL JLP_VM_READIMAG(PNTR_1,NX_1,NY_1,IN_NAME,IN_COMMENTS)

        WRITE(6,23)
23      FORMAT(' Enter the limits of the subwindow: (IX1,IY1,IX2,IY2)',/,
     1         ' (Lower left corner IX1,IY1 coordinates',
     1         ' and Upper Right corner IX2,IY2 coordinates)')
        READ(5,*) IX1,IY1,IX2,IY2

21      WRITE(6,24)
24      FORMAT(' Enter the projection angle ',
     1         '(angle of the spectrum lines relative to OX): ',/,
     1         '(between -90. and +90.; 0.0 if horizontal, 90.0 if vertical): ')
        READ(5,*) P_ANGLE
        IF(P_ANGLE.LT.-90.OR.P_ANGLE.GT.90.)GOTO 21

C Getting virtual memory: output image is NX_2 by NY_2 pixels
        NX_2=1+(IX2-IX1+1)*COS(P_ANGLE*3.14159/180.)
     1        +(IY2-IY1+1)*SIN(P_ANGLE*3.14159/180.)
C Two lines: the first is the spectral flux, the second will be the wavelengths
        NY_2=2

        ISIZE=NX_2*NY_2*4
	CALL JLP_GETVM(PNTR_2,ISIZE)

	CALL XTRACT(MADRID(PNTR_1),MADRID(PNTR_2),NX_1,NY_1,NX_2,
     1              IX1,IX2,IY1,IY2,P_ANGLE,ISTAT)

C**** Output file :
        IF(ISTAT.EQ.0)THEN
          WRITE(6,*) 'Output file: '
          READ(5,10) OUT_NAME
	  WRITE(OUT_COMMENTS,1003) IN_NAME(1:12),IX1,IY1,IX2,IY2,P_ANGLE
1003      FORMAT('From ',A,'X1,Y1,X2,Y2,ANGLE:',4(1X,I5),1X,1PG10.3)
 	  CALL JLP_WRITEIMAG(MADRID(PNTR_2),NX_2,NY_2,NX_2,
     1             OUT_NAME,OUT_COMMENTS)
        ENDIF

	CALL JLP_END
	STOP
	END
C**************************************************************
C Subroutine XTRACT
C
C**************************************************************
	SUBROUTINE XTRACT(IMAGE1,IMAGE2,NX_1,NY_1,NX_2,
     1              IX1,IX2,IY1,IY2,P_ANGLE,ISTAT)
	REAL*4 IMAGE1(NX_1,*),IMAGE2(NX_2,*)
        REAL*4 P_ANGLE,STEP2,X2,ALPHA,SC_PRO,COSP,SINP
        INTEGER*4 IX1,IX2,IY1,IY2,ISTAT,II,JJ,I2,J2

C Status:
        ISTAT=0

C Initialization of IMAGE2:
	DO I=1,NX_2
          IMAGE2(I,1)=0
        END DO

C Scalar product with vector (cos(P_ANGLE), sin(P_ANGLE))
        COSP=COS(P_ANGLE*3.14159/180.)
        SINP=SIN(P_ANGLE*3.14159/180.)
	DO J=IY1,IY2
	  DO I=IX1,IX2
            SC_PRO=FLOAT(I-IX1+1)*COSP + FLOAT(J-IY1+1)*SINP
            I2=INT(SC_PRO)
            ALPHA=SC_PRO-I2
	    IMAGE2(I2,1)=IMAGE2(I2,1)+(1.-ALPHA)*IMAGE1(I,J)
	    IMAGE2(I2+1,1)=IMAGE2(I2+1,1)+ALPHA*IMAGE1(I,J)
          END DO
        END DO

	DO I=1,NX_2
C Since the wavelengths are not known at this stage,
C we fill the second line with pixel values:
C JLP99: Inversion for speckle data:
          IMAGE2(I,2)=FLOAT(NX_2-I+1)
        END DO

	RETURN
	END

ext: for


text:
      subroutine strippath(fnam2,fnam,path,iflag)

      character*80 fnam2
      character*13 fnam
      character*80 path
      integer iflag

      integer i

      iflag=0

      do 10 i=1,80
        if(fnam2(i:i).eq.'/')then
          iflag=i
        endif
10    continue

      fnam=fnam2(iflag+1:iflag+14)

      if(iflag.gt.0)then
        path=fnam2(1:iflag)
      endif

      return
      end


MODULE TEST
CONTAINS
PURE FUNCTION s2a_3(s1,s2,s3) RESULT(a)
    CHARACTER(LEN=*), INTENT(IN)             :: s1, s2, s3
    CHARACTER(LEN=4), DIMENSION(3)        :: a

  a(1)=s1; a(2)=s2; a(3)=s3
END FUNCTION
END MODULE

USE TEST
character(len=12) :: line
write(line,'(3A4)') s2a_3("a","bb","ccc")
IF (line.NE."a   bb  ccc") CALL ABORT()
END


text:
C    %W% %G%
      subroutine znosol
c
c     THIS SUBROUTINE CALCULATES THE EQUIVALENT ADMITTANCE
c     FOR THE ZINC OXIDE ARRESTOR MODEL AND THE VARIABLE
c     ADMITTANCE ELEMENTS
c
      include 'tspinc/params.inc'
      include 'tspinc/znox.inc'
      include 'tspinc/tmptbl.inc'
      include 'tspinc/lowzno.inc'
      include 'tspinc/znox2.inc'
      include 'tspinc/vy1.inc'
      include 'tspinc/param.inc'
      include 'tspinc/vrgov.inc'
      include 'tspinc/prt.inc'
      include 'tspinc/matrow.inc'
      include 'tspinc/bname.inc'
      include 'tspinc/buskv.inc'
      include 'tspinc/spare1.inc'

      common /znotmp/ newtby, nznobg
      character*8 ibus1,jbus1

      data pi/3.1415927/

      do 1000, itrr = 1,iznmax
      lowzno = nmx + 1
c
c     STORE EXISTING TOTAL LINE EQUIVALENT IN TEMPORARY
c     VARIABLES
c
      gij = zngij(itrr)
      bij = znbij(itrr)
      gii = zngii(itrr)
      bii = znbii(itrr)
      gjj = zngjj(itrr)
      bjj = znbjj(itrr)
c
c     IF A VARIABLE ADMITTANCE ELEMENT EXISTS ON THIS LINE, THEN
c     COMBINE TOTAL VARIABLE ADMITTANCE WITH EXISTING ZNO ADMITTANCE
c
      if (vygtot(itrr) .ne. 0.0 .or. vybtot(itrr) .ne. 0.0) then
          recp1=1./(vygtot(itrr)*vygtot(itrr)
     1             +vybtot(itrr)*vybtot(itrr))
      else
          recp1 = 0.0
      endif
      if (zngc(itrr) .ne. 0.0 .or. znbc(itrr) .ne. 0.0) then
         recp2 = 1./(zngc(itrr)*zngc(itrr)
     1              + znbc(itrr)*znbc(itrr))
      else
         recp2 = 0.0
      endif
      rtot = vygtot(itrr)*recp1 + zngc(itrr)*recp2
      xtot = -(vybtot(itrr)*recp1 + znbc(itrr)*recp2)
      recp1 = 1./(rtot*rtot + xtot*xtot)
      zng12(itrr) = rtot*recp1
      znb12(itrr) =-xtot*recp1
c
c     IZNOSW = 3 MEANS THERE IS NO ZNO CAPACITOR
c
      if (iznosw(itrr) .eq. 3 .and. idsw .eq. 7) go to 1000
      if (iznosw(itrr) .eq. 3) go to 700
      if (idsw .ne. 7 .and. lppwr .eq. 0) then
c
c       IF IZNBYP = 2, THIS CAPACITOR HAS BEEN BYPASSED
c
        if (iznbyp(itrr) .ne. 2) go to 100
c
c       TEMPORARY BYPASS ON ZNO REINSERT LOGIC
c
        if (0 .eq. 0) go to 710

        if ( iznosw(itrr) .eq. 1 ) go to 950
        go to 710

c       -  Next 106 lines are bypassed
c
c       REINSRT CAP. IF CURRENT .LT. .05*IPROT
c
   40   if (ctro(itrr) .gt. 0.05 *cznpro(itrr)) then
           if ( iznosw(itrr) .eq. 1 ) go to 950
           go to 710
        endif
        ibus1 = bname(iznbus(itrr))
        jbus1 = bname(jznbus(itrr))
        write (outbuf,50 ) ibus1, jbus1, tsim
   50   format ( '0', 5x, ' CAP  REINSERTED BETWEEN BUS ',a8,
     &           ' AND ', a8, ' AT ', f8.2, ' CYCLES ' )
        call prtout(1)
        izop(itrr) = 0
        znbc(itrr) = -1./znxco(itrr)
        zngc(itrr) = 0.0
        iznbyp(itrr) = 0
        go to 620
c
c       IF ZNO STATE CHANGED OVER TWO CONSEC. TIMESTEP, PRINT INFO.FOR
c       LAST STATE
c
  100   if ( iznbyp(itrr) .ne. izop(itrr)) then
           izop(itrr) = iznbyp(itrr)
           ibus1 = bname(iznbus(itrr))
           jbus1 = bname(jznbus(itrr))
           if ( iznbyp(itrr) .eq. 0 ) then
              write (outbuf, 60 ) ibus1, jbus1, tsim
              call prtout(1)
              znengt(itrr) = 0.0
              go to 115
            endif
            if (iznbyp(itrr) .eq. 1 ) write (outbuf, 65 ) ibus1, jbus1,
     &                                tsim
  60        format('0',5x,'ZINC OXIDE ARRESTOR BETWEEN BUS ',a8,
     1         ' AND ',a8,' STOPPED CONDUCTING AT ',f8.2,' CYCLES ')
  65        format('0',5x,'ZINC OXIDE ARRESTOR BETWEEN BUS ',a8,
     1         ' AND ',a8,' CONDUCTING AT ',f8.2,' CYCLES ')
            call prtout(1)
        endif
c
c       INITIALIZE ACCUM. ENERGY IF ZNO IS NONCONDUCTING AT THE
c       LAST TIME T
C
c       CALCULATE ACCUMULATED ENERGY
c
        if ( iznbyp(itrr) .eq. 0 ) znengt(itrr) = 0.0
        znengt(itrr) = znengt(itrr) + zneng(itrr)
        if (znengt(itrr) .lt. znjou(itrr)) go to 115
        deltim = tsim - zntim(itrr)
        if ( deltim .gt. enrcyc(itrr) ) go to 115
        iznbyp(itrr) = 2
        ibus1 = bname(iznbus(itrr))
        jbus1 = bname(jznbus(itrr))
        write(outbuf, 110) ibus1, jbus1, tsim
 110    format('0',5x,'ZINC OXIDE ARRESTOR BETWEEN BUS ',a8,
     1         ' AND ',a8,' BYPASSED AT ',f8.2,' CYCLES ')
        call prtout(1)
c
c       INITIALIZE ENERGY AND TIME TABLES
c
        zntim(itrr) = endt
        zneng(itrr) = 0.0
        znengt(itrr) = 0.0
        zngc(itrr) = 0.0
        znbc(itrr) = -10000.0
        go to 620
      endif
c
c     CALL ZNMAT TO OBTAIN THE VOLTAGES ACROSS THE VARIABLE ADMITTANCE
c     SECTION--EG(1)+JFG(1) AND EG(2) + JFG(2)
c
  115 call znmat(itrr)
c
c     CALCULATE CURRENT THROUGH ZNO PROTECTED CAPACITOR
c
      zng = zng12(itrr) + zng11l(itrr)
      znb = znb12(itrr) + znb11l(itrr)
      c1r=eg(1)*zng - fg(1)*znb
      c1i=eg(1)*znb + fg(1)*zng
      c2r=eg(2)*zng12(itrr) - fg(2)*znb12(itrr)
      c2i=eg(2)*znb12(itrr) + fg(2)*zng12(itrr)
      ctr=c1r-c2r
      cti=c1i-c2i
      cmag = sqrt(ctr*ctr + cti*cti)
c
c     IF CAP IS BYPASSED STORE CURRENT AND GO TO NEXT ZNO
c
      if (iznbyp(itrr) .eq. 2 ) then
         ctro(itrr) = cmag
         go to 950
      endif
c
c       COMPARE ACTUAL CURRENT WITH PROTECTION CURRENT
c
      if (cmag .lt. .98*cznpro(itrr)) then
         zntim(itrr) = endt
         zneng(itrr) = 0.0
c
c        IF ZNO WAS CONDUCTING PREVIOUSLY, RESET ADMITTANCE TO NOMINAL
c
         if (iznbyp(itrr) .eq. 1) then
            znbc(itrr) = -1./znxco(itrr)
            zngc(itrr) = 0.0
            iznbyp(itrr) = 0
            go to 620
         endif
         iznbyp(itrr) = 0
c
c        IF ONLY ZNO AND IS NONCONDUCTING DURING PREVIOUS TIMESTEP,
c        THEN BYPASS YMATRIX UPDATE BELOW
c
         if ( iznosw(itrr) .eq. 1 ) go to 950
         go to 710
      endif
c
c     ZNO IS CONDUCTING
c
      ctr = cmag/cznpro(itrr)
c
c     IF IZNBYP = 0 THE ZNO WAS NOT CONDUCTING PREVIOUSLY
c
      if (iznbyp(itrr) .eq. 0) then
         iznbyp(itrr) = 1
         zntim(itrr) = tsim
c
c        CALCULATE NEW XC AND RC USING NEWTON'S METHOD
c
c        OPTION TO BYPASS NEWTONS APPROX. SOLN.*****

         if (newtby .eq. 2 ) go to 560
         intzn = 0
         xs = 1.0/cmag - xtot
 500     intzn = intzn +1
         if (intzn .gt. 10) then
            write(errbuf(1),550)cmag,delctr
 550        format('0',5x,'ITERATION LIMIT REACHED IN ZNO CAPACITOR ',
     1      'LOGIC CMAG AND DELCTR: ',f7.4,f9.6)
            call prterr('E',1)
            go to 1000
         endif
         xc  = 0.1010 - 0.005749*ctr + 2.088* exp(-0.8566*ctr)
         xcp = -0.005749 - 2.088 * 0.8566 * exp (-0.8566* ctr)
         fpofc = xs + znxco(itrr)*(xc + xcp*ctr)
         fofc = -1./cznpro(itrr) + ctr*(xc*znxco(itrr) + xs)
         ctr1 = ctr - fofc/fpofc
         delctr = abs(ctr1 - ctr)
         if (delctr .gt. .01) then
            ctr = ctr1
            go to 500
         endif
      endif
  560 xcpu = 0.1010 - 0.005749*ctr + 2.088* exp(-0.8566*ctr)
      rcpu = 0.0745  + 0.49*exp (-0.243*ctr) - 35.0*exp (-5.0*ctr) -
     1   0.60*exp (-1.4*ctr)
      xc =  znxco(itrr)*xcpu
      rc = -znxco(itrr)*rcpu
      recip = 1./(rc*rc + xc*xc)
      gc = rc*recip
      bc = -xc*recip
      zngc(itrr) = gc
      znbc(itrr) = bc
c
c     STORE CURRENT ESTIMATE, CTR
c
      ctro(itrr) = ctr
c
c     TALLY ENERGY OVER THIS TIMESTEP ******
c
      if (idsw .ne. 7 ) zneng(itrr) = edt *cmag*cmag*rc
c
c     COMBINE TOTAL VARIABLE ADMITTANCE WITH NEW ZNO ADMITTANCE
c
 620  if (vygtot(itrr).ne. 0.0 .or. vybtot(itrr).ne. 0.0) then
          recp1=1./(vygtot(itrr)*vygtot(itrr)
     1           +vybtot(itrr)*vybtot(itrr))
      else
          recp1 = 0.0
      endif
      recp2 = 1./(zngc(itrr)*zngc(itrr)+ znbc(itrr)*znbc(itrr))
      rtot = vygtot(itrr)*recp1 + zngc(itrr)*recp2
      xtot = -(vybtot(itrr)*recp1 + znbc(itrr)*recp2)
      recp1 = 1./(rtot*rtot + xtot*xtot)
      zng12(itrr) = rtot*recp1
      znb12(itrr) =-xtot*recp1

c
c     CALL ZNMAT TO CALCULATE NEW LINE EQUIVALENT
c
  700 call znmat(itrr)
c
c     REMOVE OLD EQUIVALENT ADMITTANCE FROM Y MATRIX
c     AND ADD NEW EQUIVALENT
c
  710 ibus = iznbus(itrr)
      jbus = jznbus(itrr)
      call getmat(ibus, ii)
      iii = ii
      do i = 4, iii, 3
        if (jbus.eq.matrow(i)) go to 920
      enddo
  920 atrow(ii-1) = atrow(ii-1) -gij - gii + zngij(itrr) + zngii(itrr)
      atrow(ii) = atrow(ii) - bij - bii + znbij(itrr) + znbii(itrr)
      atrow(i+1) = atrow(i+1) + gij  - zngij(itrr)
      atrow(i+2) = atrow(i+2) + bij  - znbij(itrr)
      call putmat(ibus, ii)
      call getmat(jbus, ii)
      iii = ii
      do i = 4, iii, 3
        if (ibus.eq.matrow(i)) go to 940
      enddo
  940 atrow(ii-1) = atrow(ii-1) - gij - gjj + zngij(itrr) + zngjj(itrr)
      atrow(ii) = atrow(ii) - bij - bjj + znbij(itrr) + znbjj(itrr)
      atrow(i+1) = atrow(i+1) + gij  - zngij(itrr)
      atrow(i+2) = atrow(i+2) + bij  - znbij(itrr)
      call putmat(jbus, ii)
c
c      GET LOWEST BUS NO. FOR YMATRIX FACTORIZING
c
       lowzno = min0 (lowzno, ibus)
       lowzno = min0 (lowzno, jbus)
  950  continue
c
c      TEMP DEBUG OUTPUT ******
c
       if ( nznobg .ne. 0 ) then
          write(outbuf, 650) itrr, iznosw(itrr), iznbyp(itrr), intzn,
     1    edt, eg(1), fg(1), eg(2), fg(2) ,zneng(itrr), znengt(itrr)
  650     format (1x, 'ZNODEBUG...', 4i5, 7e13.6 )
          call prtout (1)
          write ( outbuf, 655 )
     1    zng, znb, zngc(itrr), znbc(itrr), zng12(itrr), znb12(itrr)
  655     format ( 5x, 6e13.6 )
          call prtout(1)
          write (outbuf, 656 ) cmag, ctr, cznpro(itrr)
  656     format (1x, 3e13.6 )
          call prtout(1)
       endif
1000  continue
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


!! This module is for string manipulation
module string_utils

use, intrinsic:: iso_c_binding, only: c_null_char
implicit none
contains

elemental function toLower(str)
!! convert uppercase characters to lowercase
!!
!! can be trivially extended to non-ASCII
character(*), intent(in) :: str
character(len(str)) :: toLower
character(*), parameter :: lower="abcdefghijklmnopqrstuvwxyz", &
                           upper="ABCDEFGHIJKLMNOPQRSTUVWXYZ"
integer :: i,j

toLower = str

do concurrent (i = 1:len(str))
  j = index(upper,str(i:i))
  if (j > 0) toLower(i:i) = lower(j:j)
end do

end function toLower

pure function strip_trailing_null(str) result(stripped)
!! strip trailing C null from strings

character(*), intent(in) :: str
character(:), allocatable :: stripped
integer :: i

i = len_trim(str)
if (str(i:i) == c_null_char) then
  stripped = str(:i-1)
else
  stripped = str
endif

end function strip_trailing_null

end module string_utils


text:
      PROGRAM CMPOU
C
C Define the error file, the Fortran unit number, the workstation type,
C and the workstation ID to be used in calls to GKS routines.
C
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)   ! NCGM
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=8,  IWKID=1)   ! X Windows
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=11, IWKID=1)   ! PDF
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=20, IWKID=1)   ! PostScript
C
      PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)
C
C Open GKS, Turn Clipping off
C
      CALL GOPKS (IERRF, ISZDM)
      CALL GOPWK (IWKID, LUNIT, IWTYPE)
      CALL GACWK (IWKID)
C
C Invoke demo driver
C
      CALL TCMPOU
C
C Deactivate and close workstation, close GKS.
C
      CALL GDAWK (IWKID)
      CALL GCLWK (IWKID)
      CALL GCLKS

      STOP
      END

      SUBROUTINE TCMPOU

      REAL PLIM1(2), PLIM2(2), PLIM3(2), PLIM4(2)

      DATA PLIM1 /30.,0./
      DATA PLIM2 /-15.,0./
      DATA PLIM3 /60.,0./
      DATA PLIM4 /30.,0./
C
C CMPOU demonstrates political boundaries in the Maps utility.
C
C Set up Maps.
C
      CALL MAPSTC ('OU - OUTLINE DATASET SELECTOR','PO')
      CALL MAPROJ ('ME',0.,0.,0.)
      CALL MAPSET ('CO',PLIM1,PLIM2,PLIM3,PLIM4)
C
C Initialize Maps.
C
      CALL MAPINT
C
C Draw a perimeter and outline all the countries.
C
      CALL MAPSTI('LA - LABEL FLAG',0)
      CALL MAPLBL
      CALL MAPLOT
C
C Advance the frame.
C
      CALL FRAME
C
C Done.
C
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


subroutine load_filename_lag_ene()
  implicit none

  character (len=200) :: name_path, filename_ch

  call execute_command_line('ls')
  write(*,*)

  if (yes_no('   Have you created the initial file init.txt? ')) then
     goto 666
  else
     write(*,*) '    Create the file init.txt with this structure:'
     write(*,*) '      line 1 path name: the full path of the folder where the light curves are (without / at the end)'
     write(*,*) '      line 2 filename of the energy bins: this is where the energy bins are '
     write(*,*)
     write(*,*)

666  continue

end subroutine load_filename_lag_ene


text:
      double complex function wzdotc(n, zx, incx, zy, incy)
      double complex zx(*), zy(*), z
      double complex zdotc
      integer n, incx, incy

      z = zdotc(n, zx, incx, zy, incy)
      wzdotc = z

      end

      double complex function wzdotu(n, zx, incx, zy, incy)
      double complex zx(*), zy(*), z, zdotu
      integer n, incx, incy

      z = zdotu(n, zx, incx, zy, incy)
      wzdotu = z

      return
      end

      complex function wcdotc(n, cx, incx, cy, incy)
      complex cx(*), cy(*), c, cdotc
      integer n, incx, incy

      c = cdotc(n, cx, incx, cy, incy)
      wcdotc = c

      return
      end

      complex function wcdotu(n, cx, incx, cy, incy)
      complex cx(*), cy(*), c, cdotu
      integer n, incx, incy

      c = cdotu(n, cx, incx, cy, incy)
      wcdotu = c

      return
      end

      real function wsdot(n, x, incx, y, incy)
      real x(*), y(*), s, sdot
      integer n, incx, incy

      s = sdot(n, x, incx, y, incy)
      wsdot = s

      return
      end

      complex function wcladiv(x, y)
      complex x, y, z
      complex cladiv

      z = cladiv(x, y)
      wcladiv = z
      return
      end

      double complex function wzladiv(x, y)
      double complex x, y, z
      double complex zladiv

      z = zladiv(x, y)
      wzladiv = z
      return
      end

max_stars_repo_name: seberg/scipy

max_issues_repo_name: seberg/scipy

max_forks_repo_name: seberg/scipy


text:
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C   FFTPACK 5.0
C   Copyright (C) 1995-2004, Scientific Computing Division,
C   University Corporation for Atmospheric Research
C   Licensed under the GNU General Public License (GPL)
C
C   Authors:  Paul N. Swarztrauber and Richard A. Valent
C
C   $Id: cmfm1b.f,v 1.2 2004/06/15 21:08:32 rodney Exp $
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE CMFM1B (LOT,JUMP,N,INC,C,CH,WA,FNF,FAC)
      COMPLEX       C(*)
      REAL       CH(*),     WA(*),     FAC(*)
C
C FFTPACK 5.0 auxiliary routine
C
      NF = FNF
      NA = 0
      L1 = 1
      IW = 1
      DO 125 K1=1,NF
         IP = FAC(K1)
         L2 = IP*L1
         IDO = N/L2
         LID = L1*IDO
         NBR = 1+NA+2*MIN(IP-2,4)
         GO TO (52,62,53,63,54,64,55,65,56,66),NBR
   52    CALL CMF2KB (LOT,IDO,L1,NA,C,JUMP,INC,CH,1,LOT,WA(IW))
         GO TO 120
   62    CALL CMF2KB (LOT,IDO,L1,NA,CH,1,LOT,C,JUMP,INC,WA(IW))
         GO TO 120
   53    CALL CMF3KB (LOT,IDO,L1,NA,C,JUMP,INC,CH,1,LOT,WA(IW))
         GO TO 120
   63    CALL CMF3KB (LOT,IDO,L1,NA,CH,1,LOT,C,JUMP,INC,WA(IW))
         GO TO 120
   54    CALL CMF4KB (LOT,IDO,L1,NA,C,JUMP,INC,CH,1,LOT,WA(IW))
         GO TO 120
   64    CALL CMF4KB (LOT,IDO,L1,NA,CH,1,LOT,C,JUMP,INC,WA(IW))
         GO TO 120
   55    CALL CMF5KB (LOT,IDO,L1,NA,C,JUMP,INC,CH,1,LOT,WA(IW))
         GO TO 120
   65    CALL CMF5KB (LOT,IDO,L1,NA,CH,1,LOT,C,JUMP,INC,WA(IW))
         GO TO 120
   56    CALL CMFGKB (LOT,IDO,IP,L1,LID,NA,C,C,JUMP,INC,CH,CH,1,
     1     LOT,WA(IW))
         GO TO 120
   66    CALL CMFGKB (LOT,IDO,IP,L1,LID,NA,CH,CH,1,LOT,C,C,
     1     JUMP,INC,WA(IW))
  120    L1 = L2
         IW = IW+(IP-1)*(IDO+IDO)
         IF(IP .LE. 5) NA = 1-NA
  125 CONTINUE
      RETURN
      END

max_stars_repo_name: rpepr/rpEPR

max_issues_repo_name: rpepr/rpEPR

max_forks_repo_name: rpepr/rpEPR


text:
module bfgs

  !--------------------------------------------------------------------!
  ! Function minimization using the BFGS quasi-Newton method.          !
  !--------------------------------------------------------------------!
  ! The code is losely based on Bernd Meyer's `relax' program.         !
  !--------------------------------------------------------------------!
  ! 2011-04-08 Alexander Urban (AU)                                    !
  !--------------------------------------------------------------------!

  implicit none

  integer,          parameter, private :: u_save   = 40
  character(len=*), parameter, private :: savefile = 'BFGS.save'

  !--------------------------------------------------------------------!

  double precision, dimension(3), public  :: bfgs_a0   = 1.0d-1
  double precision,               public  :: bfgs_up   = 5.0d-5
  logical,                        public  :: bfgs_conv = .false.

  !--------------------------------------------------------------------!

  double precision, dimension(:),   allocatable, private :: bfgs_X, bfgs_X_old
  double precision, dimension(:),   allocatable, private :: bfgs_F, bfgs_F_old
  double precision, dimension(:,:), allocatable, private :: bfgs_A
  integer,                                       private :: bfgs_DoF
  logical,                                       private :: bfgs_restart

  logical, private :: isInit = .false.

contains

  subroutine bfgs_init(DoF, alpha, update, file)

    implicit none

    integer,                    intent(in) :: DoF
    double precision, optional, intent(in) :: alpha
    double precision, optional, intent(in) :: update
    logical,          optional, intent(in) :: file

    if (isInit) then
       write(0,*) "Warning: redundant initialization (module `bfgs')."
       return
    end if

    if (present(alpha)) then
       bfgs_a0(:) = alpha
    end if

    if (present(update)) then
       bfgs_up = alpha
    end if

    if ((.not. present(file)) .or. (.not. file)) then

       bfgs_DoF     = DoF
       bfgs_restart = .false.

       allocate(bfgs_X(DoF), bfgs_X_old(DoF), &
                bfgs_F(DoF), bfgs_F_old(DoF), &
                bfgs_A(DoF,DoF))

       isInit = .true.

    end if

  end subroutine bfgs_init

  !--------------------------------------------------------------------!

  subroutine bfgs_final()

    implicit none

    if (isInit) then
       deallocate(bfgs_X, bfgs_X_old, bfgs_F, bfgs_F_old, bfgs_A)
       isInit = .false.
    end if

  end subroutine bfgs_final

  !--------------------------------------------------------------------!
  !                    restart information handling                    !
  !--------------------------------------------------------------------!

  subroutine bfgs_save(n, X, F, A)

    implicit none

    integer,                          intent(in) :: n
    double precision, dimension(n),   intent(in) :: X
    double precision, dimension(n),   intent(in) :: F
    double precision, dimension(n,n), intent(in) :: A

    if (isInit)  then
       bfgs_DoF     = n
       bfgs_X(:)    = X(:)
       bfgs_F(:)    = F(:)
       bfgs_A(:,:)  = A(:,:)
       bfgs_restart = .true.
    else
       open(u_save, file=savefile, status='replace', action='write', &
                    form='unformatted')
       write(u_save) n
       write(u_save) X(:)
       write(u_save) F(:)
       write(u_save) A(:,:)
       close(u_save)
    end if

  end subroutine bfgs_save

  !--------------------------------------------------------------------!

  subroutine bfgs_restore(n, X, F, A, restart)

    implicit none

    integer,                          intent(in)  :: n
    double precision, dimension(n),   intent(out) :: X
    double precision, dimension(n),   intent(out) :: F
    double precision, dimension(n,n), intent(out) :: A
    logical,                          intent(out) :: restart

    logical :: fexists
    integer :: n1, i

    locale : if (isInit) then
       ! load from memory:
       if (n /= bfgs_DoF) then
          write(0,*) "Error: incompatible restart information" &
                  // " in `bfgs_restore()'."
          stop
       else
          if (bfgs_restart) then
             X(:)    = bfgs_X(:)
             F(:)    = bfgs_F(:)
             A(:,:)  = bfgs_A(:,:)
             restart = .true.
          else
             X(:) = 0.0d0
             F(:) = 0.0d0
             ! Initialize inverse Jacobian matrix:
             ! FIXME: generalize this to arbitrary dimensions.
             A(1:n,1:n) = 0.D0
             do i = 1, n
                A(i,i) = bfgs_a0(mod(i-1,3)+1)
             end do
             restart = .false.
          end if
       end if
    else
       inquire(file=savefile, exist=fexists)
       if (.not. fexists) then
          X(:) = 0.0d0
          F(:) = 0.0d0
          ! Initialize inverse Jacobian matrix:
          ! FIXME: generalize this to arbitrary dimensions.
          A(1:n,1:n) = 0.D0
          do i = 1, n
             A(i,i) = bfgs_a0(mod(i-1,3)+1)
          end do
          restart = .false.
       else
          ! load from file:
          open(u_save, file=savefile, status='old', action='read', &
                       form='unformatted')
          read(u_save) n1
          if (n /= n1) then
             write(0,*) "Error: incompatible restart information" &
                     // " in `bfgs_restore()'."
             stop
          else
             read(u_save) X(:)
             read(u_save) F(:)
             read(u_save) A(:,:)
          end if
          close(u_save)
          restart = .true.
       end if
    end if locale

  end subroutine bfgs_restore

  !--------------------------------------------------------------------!
  !                implementation of the BFGS algorithm                !
  !--------------------------------------------------------------------!

  subroutine bfgs_optimize(n, X, F, tol, F_max, F_rms, conv)

    implicit none

    integer,                        intent(in)    :: n
    double precision, dimension(n), intent(inout) :: X
    double precision, dimension(n), intent(in)    :: F
    double precision,               intent(in)    :: tol
    double precision,               intent(out)   :: F_max
    double precision,               intent(out)   :: F_rms
    logical,                        intent(out)   :: conv

    double precision, dimension(n)   :: X_old, F_old
    double precision, dimension(n,n) :: A
    logical                          :: restart

    call bfgs_fstatus(n, F, F_max, F_rms)
    if (abs(F_max) < tol) then
       conv = .true.
       return
    else
       conv = .false.
    end if

    call bfgs_restore(n, X_old, F_old, A, restart)

    if (restart .and. (F_rms > bfgs_up)) then
       call bfgs_jacobian(n, X, X_old, F, F_old, A)
    end if

    call bfgs_save(n, X, F, A)

    call bfgs_parameters(n, X, F, A)

  end subroutine bfgs_optimize

  !--------------------------------------------------------------------!
  !                   update inverse Jacobian matrix                   !
  !--------------------------------------------------------------------!

  subroutine bfgs_jacobian(n, X, X_old, F, F_old, A)

    implicit none

    integer,                          intent(in)    :: n
    double precision, dimension(n),   intent(in)    :: X, X_old
    double precision, dimension(n),   intent(in)    :: F, F_old
    double precision, dimension(n,n), intent(inout) :: A

    double precision              :: fac, F_norm, F_scal
    double precision, dimension(n):: av, va
    double precision, dimension(n):: dX
    double precision, dimension(n):: dF
    integer                       :: i, j

    do i = 1, n
       dX(i) = X(i) - X_old(i)
       dF(i) = F_old(i) - F(i)
    end do

    ! av(:) = matmul(A,dF) ?
    av(1:n) = 0.0d0
    do j = 1, n
       do i = 1, n
          av(i) = av(i) + A(i,j)*dF(j)
       end do
    end do

    F_scal = 0.0d0
    F_norm = 0.0d0
    do i = 1, n
       F_scal = F_scal + dX(i)*dF(i)
       F_norm = F_norm + dF(i)*av(i)
    end do
    F_scal = 1.0d0/F_scal
    fac    = 1.0d0/F_norm

    do i = 1, n
       va(i) = F_scal*dX(i) - fac*av(i)
    end do

    do j = 1, n
       do i = 1, n
          A(i,j) = A(i,j) &
               + F_scal*dX(i)*dX(j) - fac*av(i)*av(j) + F_norm*va(i)*va(j)
       end do
    end do

  end subroutine bfgs_jacobian

  !--------------------------------------------------------------------!
  !               update parameters / optimization step                !
  !--------------------------------------------------------------------!

  subroutine bfgs_parameters(n, X, F, A)

    implicit none

    integer,                          intent(in)    :: n
    double precision, dimension(n),   intent(inout) :: X
    double precision, dimension(n),   intent(in)    :: F
    double precision, dimension(n,n), intent(in)    :: A

    integer          :: i, j
    double precision :: shft

    do i = 1, n
       shft = 0.0d0
       do j = 1, n
          shft = shft + A(i,j)*F(j)
       end do
       X(i) = X(i) + shft
    end do

  end subroutine bfgs_parameters

  !====================================================================!
  !                                                                    !
  !                        auxiliary routines                          !
  !                                                                    !
  !====================================================================!

  !--------------------------------------------------------------------!
  !        maximum and root mean squared values of derivatives F       !
  !--------------------------------------------------------------------!

  subroutine bfgs_fstatus(n, F, F_max, F_rms)

    implicit none

    integer,                        intent(in)  :: n
    double precision, dimension(n), intent(in)  :: F
    double precision,               intent(out) :: F_max
    double precision,               intent(out) :: F_rms

    integer          :: i

    F_max = 0.0d0
    F_rms = 0.0d0
    do i = 1, n
       if (abs(F(i)) > abs(F_max)) F_max = F(i)
       F_rms = F_rms + F(i)*F(i)
    end do
    F_rms = sqrt(F_rms)

  end subroutine bfgs_fstatus

end module bfgs

max_stars_repo_name: alexurba/cftb

max_issues_repo_name: alexurba/cftb

max_forks_repo_name: alexurba/cftb


text:
subroutine amrex_probinit (init,name,namlen,problo,probhi) bind(c)

  use probdata_module
  use amrex_error_module
  use amrex_fort_module, only : rt => amrex_real
  implicit none

  integer :: init, namlen
  integer :: name(namlen)
  real(rt)         :: problo(2), probhi(2)

  integer untin,i

  namelist /fortin/ frac, &
       rho_1, rho_2, p0_base

  ! Build "probin" filename -- the name of file containing fortin namelist.
  integer, parameter :: maxlen = 256
  character probin*(maxlen)

  if (namlen .gt. maxlen) then
     call amrex_error('probin file name too long')
  end if

  do i = 1, namlen
     probin(i:i) = char(name(i))
  end do

  ! set namelist defaults here
  frac = 0.5e0_rt
  rho_1 = 1.0e0_rt
  rho_2 = 2.0e0_rt
  p0_base = 5.0e0_rt

  ! Read namelists
  untin = 9
  open(untin,file=probin(1:namlen),form='formatted',status='old')
  read(untin,fortin)
  close(unit=untin)

  ! set local variable defaults
  split(1) = frac*(problo(1)+probhi(1))
  split(2) = frac*(problo(2)+probhi(2))

  L_x = probhi(1) - problo(1)

end subroutine amrex_probinit

! ::: -----------------------------------------------------------
! ::: This routine is called at problem setup time and is used
! ::: to initialize data on each grid.
! :::
! ::: NOTE:  all arrays have one cell of ghost zones surrounding
! :::        the grid interior.  Values in these cells need not
! :::        be set here.
! :::
! ::: INPUTS/OUTPUTS:
! :::
! ::: level     => amr level of grid
! ::: time      => time at which to init data
! ::: lo,hi     => index limits of grid interior (cell centered)
! ::: nstate    => number of state components.  You should know
! :::		   this already!
! ::: state     <=  Scalar array
! ::: delta     => cell size
! ::: xlo,xhi   => physical locations of lower left and upper
! :::              right hand corner of grid.  (does not include
! :::		   ghost region).
! ::: -----------------------------------------------------------
subroutine ca_initdata(level,time,lo,hi,nscal, &
                       state,state_l1,state_l2,state_h1,state_h2, &
                       delta,xlo,xhi)

  use probdata_module
  use meth_params_module, only : NVAR, URHO, UMX, UMY, &
       UEDEN, UEINT, UFS, UTEMP, small_temp
  use amrex_constants_module, only: ZERO, HALF, M_PI
  use actual_eos_module, only : gamma_const

  use amrex_fort_module, only : rt => amrex_real
  implicit none

  integer :: level, nscal
  integer :: lo(2), hi(2)
  integer :: state_l1,state_l2,state_h1,state_h2
  real(rt)         :: xlo(2), xhi(2), time, delta(2)
  real(rt)         :: state(state_l1:state_h1,state_l2:state_h2,NVAR)

  integer :: i,j
  real(rt)         :: x,y,pres,presmid,pertheight

  presmid  = p0_base - rho_1*split(2)

  state(:,:,UMX)   = ZERO
  state(:,:,UMY)   = ZERO
  state(:,:,UTEMP) = small_temp

  do j = lo(2), hi(2)
     y = (j+HALF)*delta(2)

     do i = lo(1), hi(1)

        if (y .lt. split(2)) then
           pres = p0_base - rho_1*y
           state(i,j,UEDEN) = pres / (gamma_const - 1.0e0_rt)
           state(i,j,UEINT) = pres / (gamma_const - 1.0e0_rt)
        else
           pres = presmid - rho_2*(y-split(2))
           state(i,j,UEDEN) = pres / (gamma_const - 1.0e0_rt)
           state(i,j,UEINT) = pres / (gamma_const - 1.0e0_rt)
        end if

     enddo
  enddo

  do j = lo(2), hi(2)
     y = (j+HALF)*delta(2)

     do i = lo(1), hi(1)
        x = (i+HALF)*delta(1)

        ! we explicitly make the perturbation symmetric here
        ! -- this prevents the RT from bending.
        pertheight = 0.01e0_rt*HALF*(cos(2.0e0_rt*M_PI*x/L_x) + &
                                  cos(2.0e0_rt*M_PI*(L_x-x)/L_x)) + 0.5e0_rt
        state(i,j,URHO) = rho_1 + ((rho_2-rho_1)/2.0e0_rt)* &
             (1+tanh((y-pertheight)/0.005e0_rt))
        state(i,j,UFS) = state(i,j,URHO)

     enddo
  enddo

end subroutine ca_initdata

max_stars_repo_name: harpolea/Castro

max_issues_repo_name: harpolea/Castro

max_forks_repo_name: harpolea/Castro


text:
FUNCTION pythag_sp(a,b)
  USE nrtype
  IMPLICIT NONE
  REAL(SP), INTENT(IN) :: a,b
  REAL(SP) :: pythag_sp
  REAL(SP) :: absa,absb
  absa=abs(a)
  absb=abs(b)
  if (absa > absb) then
     pythag_sp=absa*sqrt(1.0_sp+(absb/absa)**2)
  else
     if (absb == 0.0) then
        pythag_sp=0.0
     else
        pythag_sp=absb*sqrt(1.0_sp+(absa/absb)**2)
     end if
  end if
END FUNCTION pythag_sp

FUNCTION pythag_dp(a,b)
  USE nrtype
  IMPLICIT NONE
  REAL(DP), INTENT(IN) :: a,b
  REAL(DP) :: pythag_dp
  REAL(DP) :: absa,absb
  absa=abs(a)
  absb=abs(b)
  if (absa > absb) then
     pythag_dp=absa*sqrt(1.0_dp+(absb/absa)**2)
  else
     if (absb == 0.0) then
        pythag_dp=0.0
     else
        pythag_dp=absb*sqrt(1.0_dp+(absa/absb)**2)
     end if
  end if
END FUNCTION pythag_dp

max_stars_repo_name: aasensio/DeepLearning

max_issues_repo_name: aasensio/DeepLearning

max_forks_repo_name: aasensio/DeepLearning


! { dg-do compile }
! { dg-compile-aux-modules "binding_label_tests_13.f03" }
module binding_label_tests_13_main
  use, intrinsic :: iso_c_binding, only: c_int
  integer(c_int) :: c3  ! { dg-error "Variable c3 from module binding_label_tests_13_main with binding label c3 at .1. uses the same global identifier as entity at .2. from module binding_label_tests_13" }
  bind(c) c3

contains
  subroutine c_sub() BIND(c, name = "C_Sub")
    use binding_label_tests_13 ! { dg-error "Variable c3 from module binding_label_tests_13_main with binding label c3 at .1. uses the same global identifier as entity at .2. from module binding_label_tests_13" }
  end subroutine c_sub
end module binding_label_tests_13_main
! { dg-final { cleanup-modules "binding_label_tests_13" } }


text:
! MIT License
!
! Copyright (c) 2010-present David A. Kopriva and other contributors: AUTHORS.md
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!
! HOHQMesh contains code that, to the best of our knowledge, has been released as
! public domain software:
! * `b3hs_hash_key_jenkins`: originally by Rich Townsend,
!    https://groups.google.com/forum/#!topic/comp.lang.fortran/RWoHZFt39ng, 2005
! * `fmin`: originally by George Elmer Forsythe, Michael A. Malcolm, Cleve B. Moler,
!    Computer Methods for Mathematical Computations, 1977
! * `spline`: originally by George Elmer Forsythe, Michael A. Malcolm, Cleve B. Moler,
!    Computer Methods for Mathematical Computations, 1977
! * `seval`: originally by George Elmer Forsythe, Michael A. Malcolm, Cleve B. Moler,
!    Computer Methods for Mathematical Computations, 1977
!
! --- End License
!
!////////////////////////////////////////////////////////////////////////
!
!      Filename.f
!      Created: 2010-08-27 12:14:31 -0400
!      By: David Kopriva
!
! /////////////////////////////////////////////////////////////////////
!
   SUBROUTINE toLower(str)
!
!  ----------------
!  From ResettaCode
!  ----------------
!
     character(*), intent(in out) :: str
     integer :: i

     do i = 1, len(str)
       select case(str(i:i))
         case("A":"Z")
           str(i:i) = achar(iachar(str(i:i))+32)
       end select
     end do
   END SUBROUTINE toLower
!
!///////////////////////////////////////////////////////////////////////
!
!     ----------------------------------------------------------------
!>    "Read" the "value" of a real number declared
!!     after an = sign in a inputLine
!     ----------------------------------------------------------------
!
      FUNCTION GetRealValue( inputLine )
         USE, INTRINSIC :: iso_fortran_env, only : stderr => ERROR_UNIT
         USE SMConstants
         IMPLICIT NONE
!
         CHARACTER(len = *) :: inputLine
         REAL( KIND=RP )    :: value, GetRealValue
         INTEGER            :: strLen, leq, k
!
         leq    = INDEX( inputLine, '=' )
         strLen = LEN_TRIM( inputLine )
         READ( inputLine( leq+1:strLen ), *, IOSTAT=k ) VALUE
         IF ( k /=0 )     THEN
            WRITE(stderr,*)  "Bad real value in input line:"
            WRITE(stderr,*)  TRIM(inputLine)
            ERROR STOP "Input file synax error"
         END IF
         GetRealValue = value
!
      END FUNCTION GetRealValue
!
!///////////////////////////////////////////////////////////////////////
!
!     ----------------------------------------------------------------
!>    "Read" the "value" of an integer number declared
!!     after an = sign in a inputLine
!     ----------------------------------------------------------------
!
      FUNCTION GetRealArray( inputLine ) RESULT(x)
         USE SMConstants
         USE, INTRINSIC :: iso_fortran_env, only : stderr => ERROR_UNIT
         IMPLICIT NONE
!
         REAL(KIND=RP), DIMENSION(3) :: x
         INTEGER                     :: k

         CHARACTER ( LEN = * ) :: inputLine
         INTEGER               :: cStart, cEnd
!
         cStart = INDEX(inputLine,'[')
         cEnd   = INDEX(inputLine, ']', .true. )
         READ( inputLine( cStart+1: cEnd-1 ), *, IOSTAT = k  ) x(1), x(2), x(3)
         IF ( k /=0 )     THEN
            WRITE(stderr,*)  "Bad real array syntax in input line:"
            WRITE(stderr,*)  TRIM(inputLine)
            WRITE(stderr,*)  "Syntax is: [real,real,real]"
            ERROR STOP "Input file synax error"
         END IF
!
      END FUNCTION GetRealArray
!
!///////////////////////////////////////////////////////////////////////
!
!     ----------------------------------------------------------------
!>    "Read" the "value" of an integer number declared
!!     after an = sign in a inputLine
!     ----------------------------------------------------------------
!
      INTEGER FUNCTION GetIntValue( inputLine )
         USE, INTRINSIC :: iso_fortran_env, only : stderr => ERROR_UNIT
         IMPLICIT NONE
!
         CHARACTER ( LEN = * ) :: inputLine
         INTEGER               :: value
         INTEGER               :: strLen, leq, k
!
         leq    = INDEX( inputLine, '=' )
         strLen = LEN_TRIM( inputLine )
         READ( inputLine( leq+1:strLen ), *, IOSTAT = k ) VALUE
         IF ( k /=0 )     THEN
            WRITE(stderr,*) "Bad integer value in input line:"
            WRITE(stderr,*) TRIM(inputLine)
            ERROR STOP "Input file synax error"
         END IF
         GetIntValue = VALUE
!
      END FUNCTION GetIntValue
!
!///////////////////////////////////////////////////////////////////////
!
!     ----------------------------------------------------------------
!>    "Read" the "value" of an integer number declared
!!     after an = sign in a inputLine
!     ----------------------------------------------------------------
!
      FUNCTION GetIntArray( inputLine ) RESULT(N)
         USE, INTRINSIC :: iso_fortran_env, only : stderr => ERROR_UNIT
         IMPLICIT NONE
!
         INTEGER, DIMENSION(3) :: N

         CHARACTER ( LEN = * ) :: inputLine
         INTEGER               :: cStart, cEnd, k
!
         cStart = INDEX(inputLine,'[')
         cEnd   = INDEX(inputLine, ']', .true. )
         READ( inputLine( cStart+1: cEnd-1 ), *, IOSTAT = k ) N(1), N(2), N(3)
         IF ( k /=0 )     THEN
            WRITE(stderr,*)  "Bad integer array value in input line:"
            WRITE(stderr,*)  TRIM(inputLine)
            WRITE(stderr,*)  "Syntax is: [integer,integer,integer]"
            ERROR STOP "Input file synax error"
         END IF
!
      END FUNCTION GetIntArray
!
!///////////////////////////////////////////////////////////////////////
!
!     ----------------------------------------------------------------
!>    Extracts the string within the quotes in an input file
!     ----------------------------------------------------------------
!
      CHARACTER( LEN=LINE_LENGTH ) FUNCTION GetStringValue( inputLine )
         USE ProgramGlobals
         IMPLICIT NONE
!
         CHARACTER ( LEN = * ) :: inputLine
         INTEGER               :: cStart, cEnd
!
         cStart = INDEX(inputLine,'"')
         cEnd   = INDEX(inputLine, '"', .true. )
         GetStringValue = inputLine( cStart+1: cEnd-1 )
!
      END FUNCTION GetStringValue
!
!///////////////////////////////////////////////////////////////////////
!
!     ----------------------------------------------------------------
!>    "Read" the "value" of an logical declared
!!     after an = sign in a inputLine
!     ----------------------------------------------------------------
!
      LOGICAL FUNCTION GetLogicalValue( inputLine )
         USE, INTRINSIC :: iso_fortran_env, only : stderr => ERROR_UNIT
         IMPLICIT NONE
!
         CHARACTER ( LEN = * ) :: inputLine
         LOGICAL               :: value
         INTEGER               :: strLen, leq, k
!
         leq    = INDEX( inputLine, '=' )
         strLen = LEN_TRIM( inputLine )
         READ( inputLine( leq+1:strLen ), *, IOSTAT = k ) value
         GetLogicalValue = value
         IF ( k /=0 )     THEN
            WRITE(stderr,*)  "Bad logical value in input line:"
            WRITE(stderr,*)  TRIM(inputLine)
            ERROR STOP "Input file synax error"
         END IF
!
      END FUNCTION GetLogicalValue
!
! /////////////////////////////////////////////////////////////////////
!
!     ----------------------------------------------------------------
!>    This subroutine returns an unused
!!    unit number. It is assumed that the unit number will be used
!!    immediately after it is assigned. It is also assumed that the
!!    compiler has access to units 1-99. The function returns "0"
!!    if it cannot find an unused unit number.
!     ----------------------------------------------------------------
!
      INTEGER FUNCTION UnusedUnit()
!
         IMPLICIT NONE
         INTEGER :: j
         LOGICAL :: unitIsOpened, unitDoesExist

         UnusedUnit = 0
         DO j = 1, 99
           INQUIRE( UNIT = j, OPENED = unitIsOpened, EXIST = unitDoesExist )
           IF ( .NOT.unitIsOpened ) EXIT
         END DO

         IF (  j <= 99 .AND. unitDoesExist )     THEN
            UnusedUnit = j
         ELSE
            UnusedUnit = 0
         END IF
!
      END FUNCTION UnusedUnit
!
!////////////////////////////////////////////////////////////////////////
!
      INTEGER FUNCTION StdInFileUnitCopy( ) RESULT(fUnit)
!
!     ----------------------------------------------------
!     Since StdIn cannot be rewound, make a temporary file
!     ----------------------------------------------------
!
         IMPLICIT NONE
         INTEGER, EXTERNAL  :: UnusedUnit
         CHARACTER(LEN=132) :: inputLine

         fUnit = UnusedUnit( )
         OPEN(UNIT=fUnit, STATUS="scratch")

         DO
           READ ( 5, FMT = '(a132)', END = 1000 ) inputLine
           IF( INDEX(inputline, "\end{FILE}") /= 0 )     EXIT
           WRITE( fUnit,  FMT = '(a132)' ) inputLine
         END DO
 1000 CONTINUE
      rewind(fUnit)
      END FUNCTION StdInFileUnitCopy
!
!////////////////////////////////////////////////////////////////////////
!
      SUBROUTINE ConvertToPath(str)
         IMPLICIT NONE
         CHARACTER(LEN=*) :: str
         INTEGER          :: lc

         lc = LEN_TRIM(str)
         IF ( str(lc:lc) == "/" )     THEN
            str = str(1:lc)
         ELSE
            str = str(1:lc) // "/"
         END IF

      END SUBROUTINE ConvertToPath

max_stars_repo_name: HigherOrderMethods/HOHQMesh

max_issues_repo_name: DavidAKopriva/HOHQMesh

max_forks_repo_name: DavidAKopriva/HOHQMesh