program usuki
	use IMSL
	implicit none
	integer, parameter::rows=64				
	integer, parameter::cols=55				
	real(8), parameter::a=5.0d-9			
	real(8), parameter::b=0.0d0				
	real(8), parameter::EF0=0.02				
	complex(8) H(rows,rows),cimag,phase,evs(rows),evec(rows,rows)
	complex(8) inter(rows,rows),LPINV(rows,rows),H2(2*rows,2*rows)
	complex(8) k(rows),CM1(rows,rows),CM2(rows,rows),T21(rows,rows)
	complex(8) devec(2*rows,2*rows),devs(2*rows),UM(rows,rows)
	complex(8) T12(rows,rows),LP(rows,rows),inter2(rows,rows)
	complex(8) P2S(rows,rows,cols),P1S(rows,rows,cols),PHI1(rows,rows)
	complex(8) PHI2(rows,rows)
	real(8) c0,c1,c2,thop,mass,bet,hb2o2m,vel(rows),EF,bl,sum,test
	real(8) test2,dens(rows,cols),pot(rows,cols)
	real(8) hbar,m0,q
	integer i,j,s,s2,mode
	open (3,file='eigenvalues.txt')
	open (4,file='density-55.txt')

c	Define the constants

	m0=9.11d-31
	hbar=1.054d-34	
	q=1.602d-19
	mass=.067d0*m0
	c0=0.0d0
	c1=1.0d0
	c2=2.0d0
	cimag=dcmplx(c0,c1)
	hb2o2m=(hbar/c2/mass)*(hbar/q)
	thop=hb2o2m/a/a
	bet=q*b*a*a/hbar

c	Determine the eigenvectors and eigenvalues of the -1 column

	H=dcmplx(c0,c0)
	pot=c0
	CM2=H
	CM1=H
	LPINV=H
	LP=H
	do i=1,rows
		H(i,i)=c2*thop
		if (i.lt.rows) then
			H(i,i+1)=thop
		endif
		if (i.gt.rows) then
			H(i,i-1)=thop
		endif
	enddo
	call devccg(rows,H,rows,evs,evec,rows)

c	The matrix evecn now contains the eigenvectors of the zero 
c	column as the column with the lowest eigenvector in the
c	rightmost column.  We are in the site representation, 
c	but evec provides the transformation into this representation 
c	from the mode representation.

c	Ramp the fermi energy

	do i=1,rows
		write (3,*) real(evs(i))
	enddo

	test=real(evs(rows))

c	Define the local potential in terms of the hopping term

	do j=21,31
		do i=1,26
			pot(i,j)=c2*c2
			pot(64+1-i,j)=c2*c2
	enddo
	enddo

c	Initialize a loop over a range of Fermi energies

	do s=16,16
		EF=EF0*(dble(s)/20.0d0)+test

c	Initialize the modes and find the new eigenvalues for the 0 column

		H2=dcmplx(c0,c0)
		do i=1,rows
			H2(i,i+rows)=c1
			H2(i+rows,i)=-c1
			H2(i+rows,i+rows)=-(EF/thop)+c2*c2
			if (i.lt.rows) then
				H2(i+rows,i+1+rows)=-c1
			endif
			if (i.gt.1) then
				H2(i+rows,i-1+rows)=-c1
			endif
		enddo
		call devccg(2*rows,H2,2*rows,devs,devec,2*rows)
		j=1
		do i=1,2*rows
			if (real(devs(i)).gt.c1) then
				evs(j)=devs(2*rows-i+1)
				LPINV(j,j)=devs(i)
				LP(j,j)=c1/devs(i)
				do s2=1,rows
					evec(s2,j)=devec(s2,i)
				enddo
				j=j+1
			else if (imag(devs(i)).gt.1.0d-05) then
				evs(j)=devs(i)
				LP(j,j)=devs(i)
				LPINV(j,j)=dconjg(devs(i))
				do s2=1,rows
					evec(s2,j)=devec(s2,i)
				enddo
			endif
		enddo

c	Determine the propagation vectors and mode velocities

		do i=1,2*rows
			if (dabs(dimag(evs(i))).lt.1.0d-10) then
				k(i)=dcmplx(c0,-log(real(evs(i)))/a)
				vel(i)=c0
			else
				k(i)=dcmplx(atan(imag(evs(i))/real(evs(i)))/a,c0)
				vel(i)=(hbar/c2/a/mass)*dsin(c2*dreal(k(i))*a)
			endif
		enddo

c	Normalize the eigenvectors

		do i=1,rows
			test2=c0
			do j=1,rows
				test2=test2+a*cdabs(evec(j,i))**2
			enddo
			test2=dsqrt(test2)
			do j=1,rows
				evec(j,i)=evec(j,i)/test2
			enddo
		enddo
		UM=evec
		inter=matmul(UM,LPINV)

c	Compute the inverse of the matrix

		call dlincg(rows,inter,rows,CM1,rows)
		CM2=matmul(UM,CM1)
		inter=matmul(evec,LP)
		CM1=evec-matmul(CM2,inter)

c	Ramp the magnetic field over 10 columns

		do j=1,10
			bl=dble(j)*bet/10.0d0
			H=dcmplx(c0,c0)
			T21=dcmplx(c0,c0)
				do i=1,rows
					phase=cdexp(-cimag*bl*dble(i))
					H(i,i)=-((EF/thop)-c2*c2)*phase
					T21(i,i)=-phase*phase
					if (i.lt.rows) then
						H(i,i+1)=-phase
					endif
					if (i.gt.1) then 
						H(i,i-1)=-phase
					endif
				enddo
			inter=matmul(T21,CM2)+H

c	Compute the inverse of the matrix

			call dlincg(rows,inter,rows,CM2,rows)

			inter=matmul(T21,CM1)

			CM1=-matmul(CM2,inter)

		enddo
		
c	Now propagate through the active region

		do j=1,cols
			H=dcmplx(c0,c0)
			T21=dcmplx(c0,c0)
			do i=1,rows
				phase=cdexp(-cimag*bet*dble(i))
				H(i,i)=-((EF/thop)-c2*c2-pot(i,j))*phase
				T21(i,i)=-phase*phase
				if (i.lt.rows) then
					H(i,i+1)=-phase
				endif
				if (i.gt.1) then
					H(i,i-1)=-phase
				endif
		enddo
		inter=matmul(T21,CM2)+H

c	Now compute the inverse of the matrix again

		call dlincg(rows,inter,rows,CM2,rows)

		inter=matmul(T21,CM1)

		CM1=-matmul(CM2,inter)		

		P2S(:,:,j)=CM2

		P1S(:,:,j)=CM1

		enddo

c	Ramp down the magnetic field

		do j=0,10	
			bl=dble(10-j)*bet/10.0d0
			H=dcmplx(c0,c0)
			T21=dcmplx(c0,c0)
			do i=1,rows
				phase=cdexp(-cimag*bl*dble(i))
				H(i,i)=-((EF/thop)-c2*c2)*phase
				T21(i,i)=-phase*phase
				if (i.lt.rows) then
					H(i,i+1)=-phase
				endif
				if (i.gt.1) then
					H(i,i-1)=-phase
				endif
			enddo
			inter=matmul(T21,CM2)+H
			call dlincg(rows,inter,rows,CM2,rows)
			inter=matmul(T21,CM1)
			CM1=-matmul(CM2,inter)
			if (j.eq.0) then
				PHI1=CM1
				PHI2=CM2
			else
				PHI1=PHI1+matmul(PHI2,CM1)
				PHI2=matmul(PHI2,CM2)
			endif
		enddo

c	Now compute the transmission

		inter=matmul(evec,LP)
		call dlincg(rows,inter,rows,T12,rows)
		inter=CM2-matmul(evec,T12)
		call dlincg(rows,inter,rows,inter2,rows)
		inter=-matmul(inter2,CM1)
		PHI1=PHI1+matmul(PHI2,inter)
		PHI2=matmul(PHI2,inter2)
		CM2=matmul(T12,inter2)
		CM1=matmul(T12,inter)
		sum=c0
		do i=1,rows
			do j=1,rows
				if (vel(j).gt.c0) then
					sum=sum+(vel(i)*cdabs(CM1(i,j))**2)/vel(j)
				endif
			enddo
		enddo
		write(3,*) 'Transmission=',sum

c	Now unfold the density

		do j=cols,1,-1
			PHI1=P1S(:,:,j)+matmul(P2S(:,:,j),PHI1)
			do i=1,rows
				sum=c0
				do mode=1,rows
					sum=sum+cdabs(PHI1(i,mode))**2
				enddo
				dens(i,j)=sum
			enddo
		enddo
		do j=1,cols
			do i=1,rows
				write(4,*) dens(i,j)
			enddo
		enddo
		do i=1,rows
			write(3,*) k(i),vel(i)
		enddo
	enddo
	end