program test_generate_transformer
	
		Implicit None

!\\\\\\\\\\\\\\\		
		Integer :: StatusOne,StatusTwo,StatusThree,StatusFour, elleMax, nSquared, elles, i
		double precision :: time, Tau,alphaPerpendicular, alphaParallel,radius, frequency, phase, amplitude, protonMass, delta
		double precision, allocatable, dimension(:,:) :: H_Matrix
		complex, allocatable, dimension(:,:) :: transformer
!///////////////

		
!\\\\\\\\\\\\\\\			
		elleMax = 10	  
	  	elles = (elleMax+1)
		
		Allocate(H_Matrix(elles,elles), Stat = statusOne) 
		Write(*,*) 'The status will be 0 for successful memory Allocation. STATUS1 = ', StatusOne
		
		call import_mathematica_data(elleMax, elles, H_Matrix)
!///////////////		
		
		
!\\\\\\\\\\\\\\\			
		time = 6.789
		Tau = 5
		alphaPerpendicular = 0.175165E1		!These are the values at "equilibrium"
		alphaParallel = 0.507765E1		!These are the values at "equilibrium"
		Radius = 2.0				!These are the values at "equilibrium"
		frequency = 0.1E15
		phase = 0.0d0
		amplitude = 50
		protonMass = 1836.15		
		
		call generate_H(H_Matrix, elles, time, Tau, alphaPerpendicular, alphaParallel, frequency, phase, amplitude, protonMass, Radius)
!///////////////	
	
	
!\\\\\\\\\\\\\\\	
		delta = 0.5
		allocate(transformer(elles,elles), Stat = StatusTwo)
		Write(*,*) 'The status will be 0 for successful memory Allocation. STATUS2 = ', StatusTwo
		call generate_transformer(delta, time, H_Matrix, transformer, elles)
		do i = 1,elles
            		write(6,*) transformer(i,1:elles)
          	end do
		Deallocate(transformer, Stat = StatusThree)
		Write(*,*) 'The status will be 0 for successful memory DEallocation. STATUS3 = ', StatusFour
!///////////////


!\\\\\\\\\\\\\\\	
		do i = 1,elles
            		write(6,*) H_Matrix(i,1:elles)
          	end do
		
		Deallocate(H_Matrix, Stat = StatusFour)
		Write(*,*) 'The status will be 0 for successful memory DEallocation. STATUS4 = ', StatusFour
!///////////////	
	
	Stop
	End Program test_generate_transformer
	
	Subroutine generate_transformer(stepSize, currentTime, H_in, out_matrix, sideLength)
	
		Implicit None
		
		Integer, intent(in) :: sideLength
		double precision, Intent(in) :: stepSize, currentTime
		double precision, dimension(sideLength, sideLength), Intent(in) :: H_in
		complex, dimension(sideLength, sideLength), Intent(out) :: out_Matrix	
		
		Integer :: row,collumn,INFO, INFO2
		double precision, dimension(sideLength, sideLength) :: Identity
		complex, dimension(sideLength, sideLength) :: Intermediate_Plus, Intermediate_Minus
		complex*16, dimension(SideLength,sideLength) :: A
		integer, dimension(sideLength) :: IPIV
		complex*16, dimension(sideLength) :: WORK
		complex*16, dimension(2*sideLength) :: WORK2
		
		do row = 1,sideLength
			do collumn = 1,sideLength
				If(row==collumn)Then
					Identity(row,collumn) = 1.0
				Else
					Identity(row,collumn) = 0.0
				End If
			End do
		End do
		
		do row = 1,sideLength
			do collumn = 1,sideLength
				Intermediate_Plus(row,collumn) = (1,0)*Identity(row,Collumn)+(0,1)*H_in(row,collumn)*(StepSize/2)	
			End do
		End do
	
		do row = 1,sideLength
			do collumn = 1,sideLength
				Intermediate_Minus(row,collumn) = (1,0)*Identity(row,Collumn)-(0,1)*H_in(row,collumn)*(StepSize/2)	
			End do
		End do
		
		call zsytrf('U',sideLength,Intermediate_Plus,sideLength,IPIV,WORK,sideLength,INFO)
		
		call zsytri('U',sideLength,Intermediate_Plus,sideLength,IPIV,WORK2,INFO2 )
		
		out_Matrix = MatMul(Intermediate_Plus,Intermediate_Minus)
	
	End subroutine generate_transformer
	
	
	Subroutine generate_H(H_Matrix_SR, elles_SR, t_SR, Tau_SR, a_Perp_SR, a_par_SR, freq_SR, phase_SR, amp_SR, Mproton_SR, R_SR)
		
		Implicit None
		
		
		Integer, Intent(in) :: elles_SR
		double precision, Intent(in) :: t_SR, Tau_SR, A_Perp_SR, a_Par_SR, freq_SR, phase_SR, amp_SR, Mproton_SR, R_SR
		double precision, dimension(elles_SR,elles_SR), intent(inout) :: H_Matrix_SR
		
		Integer :: row,collumn
		double precision :: EveryElementCoefficient
		double precision, dimension(elles_SR, elles_SR) :: Identity
		
		EveryElementCoefficient =(amp_SR*Exp(-((t_SR/Tau_SR)**2))*Cos(freq_SR*t_SR+phase_SR)*(A_Perp_SR - a_Par_SR))/2
		
		H_Matrix_SR = H_Matrix_SR*EveryElementCoefficient
		
		do row = 1,elles_SR
			do collumn = 1,elles_SR
				If(row==collumn)Then
					Identity(row,collumn) = 1*(row*(row+1))*(1/Mproton_SR)*(1/(R_SR**2))
				Else
					Identity(row,collumn) = 0.0d0
				End If
			End do
		End do
		
		H_Matrix_SR = H_Matrix_SR+Identity
		
		Return
	
	End Subroutine generate_H


	
	Subroutine import_mathematica_data(elleMaxRoutine, ellesRoutine, importedMatrixRoutine)
	
		Implicit None
		
		Integer, Intent(in) :: elleMaxRoutine, ellesRoutine
		double precision, dimension(ellesRoutine,ellesRoutine), intent(out) :: importedMatrixRoutine
		
		integer :: row, collumn
		
		Open(7, File = '/lhome/theodore/linux/Hamiltonia/H_L_M_10_0.dat', Status = 'old', Action = 'Read')
		Read(7,*) ((importedMatrixRoutine(row,collumn), row = 1,ellesRoutine), collumn = 1,ellesRoutine)
		
		return
		
	End Subroutine import_mathematica_data