program hello
    use iso_fortran_env, only: int16, int32, real64
    use ifport
    implicit none

    integer, parameter :: RK=real64, Ishort=int16, Ilong=int32

    !Units: GHz, ns, J, gauss, mm
    integer(Ilong) :: nFreqs, nPhis, nPhis_overall, expNum, parLoopInd, folStat
    integer(Ishort) :: init_i, nStates=9
!    integer(Ishort) :: init_i, nStates=4
    integer(Ilong) :: i, j, f0Ind
    real(RK), dimension(:), allocatable :: phis, phis_loop, phis_overall, freqs_loop, loopArr, chirpScales
    real(RK), parameter :: pi=4.0_RK*ATAN(1.0_RK)
    real(RK), parameter :: h_planck=6.626070040E-25_RK, hbar=h_planck/2.0_RK/pi !CODATA, J/GHz
    real(RK), parameter :: muB=9.274009994E-28_RK !CODATA, J/gauss
    real(RK), parameter :: freq01=29.616951_RK, freq12=2.29117659_RK, EJ1=0.0_RK, EJ0=h_planck*freq01+EJ1, EJ2 = EJ1-h_planck*freq12
    real(RK), parameter :: gS=2.00231930436182_RK, gL=0.9998629_RK !CODATA
!    real(RK), parameter :: gS=2.00223886745683_RK, gL=0.9998736_RK !Drake (1994) Corrections
    real(RK), parameter :: chi=muB*(gS+gL), xi=muB*(gS-gL), saveStep = 0.1_RK
    real(RK) :: t_max, h, t, tT, tD, freqRF, amp=0.0_RK, phi, phi_chirp, PSlope, PwrImbalance, timeShiftScale
    real(RK) :: loopVar, centreT, Bz0, nPhis_overall_float, theta_MW, powerLvl, ampScale, chirpScale, chirpStretch
    complex(RK), parameter :: ci=(0.0_RK,1.0_RK), c0=(0.0_RK,0.0_RK), c1=(1.0_RK,0.0_RK)
    complex(RK), dimension(:,:), allocatable :: psi, psi_new, k, k1, k2, k3, k4
    real(RK), dimension(:,:), allocatable :: B, B_Step, B_halfStep, B_background
    complex(RK), dimension(:), allocatable :: psi_vec
    character(len=*), parameter :: resultsStem='I:\Fortran code\30GHz 9 state no read codeblocks\results'
    character(len=*), parameter :: WriteFmtStr='(E25.17,A,E25.17)'
    character(len=*), parameter :: WriteFmtStr2='(E25.17,A,E25.17E3,A,E25.17)'
    character(len=*), parameter :: WriteFmtStr3='(E25.17,A,E25.17,A,E25.17,A,E25.17,A,E25.17)'
    character(len=*), parameter :: WriteFmtStr4='(E25.17,SP,E25.17,"j",A,E25.17,SP,E25.17,"j",A,E25.17,SP,E25.17,"j")'
    character(len=100) :: jstr, loopStr, loopSizeStr
    character(len=1), dimension(2) :: expTypes=['A','B']
!    character(len=1), dimension(1) :: expTypes=['A']
    character(len=260) :: resultsDir, fullSaveDir
    character(len=1) :: expType

    ! AmpToBz = 2.0_RK/xi!*SQRT(95.0_RK/35.0_RK)

!    h = 5.0E-4_RK !Big step size (TOO BIG! ~20kHz shift)
!    h = 1.0E-4_RK !Small step size
    h = 5.0E-5_RK !Smaller step size
    tT = 300.0_RK !Experiment period, ns
    tD = 100.0_RK !Pulse duration, ns
    PSlope = 0*1.0_RK !0.2 Linear slope in power, percent per MHz
    PwrImbalance = 0*10.0_RK !Percent difference in power between pulses
    Bz0 = 4.762_RK !DC B field
    theta_MW = 0.0_RK !Microwave polarisation angle, degrees
    chirpStretch = 1.0_RK !Larger = wider chirp
    chirpScale = 1.0_RK !Scale size of all chirps
    chirpScales = [0, 0, 0, 0]/6.0_RK !Scale individual chirps
    timeShiftScale = 0.1_RK !Shift pulse separation by this amount, ns/MHz

!    loopArr = [0:20]*0.1_RK !Bz0 fields
!    loopArr = [0:10]*0.2_RK !Bz0 fields
!    loopArr = [1,4,7,10]*10.0_RK !Powers (%)
    loopArr = [10,55,100]*1.0_RK !Powers (%)
!    loopArr = [100]*1.0_RK !Powers (%)
    phis_overall = [0:2]*2.0_RK*pi/3.0_RK !For averaging over total phase

    allocate(psi_vec(nStates))
!    psi_vec = SQRT([c0, 0.1_RK*c1, 0.8_RK*c1, 0.1_RK*c1, c0, c0, c0, c0, c0]) !Start with mixed population
    psi_vec = SQRT([c0, c0, c1, c0, c0, c0, c0, c0, c0]) !Start with perfect population
!    psi_vec = SQRT([c0, c0, c1, c0])

!    freqs_loop = [-21:21]/3.0_RK
!    freqs_loop = [-12:12]/2.0_RK
!    freqs_loop = [-25:25]*1.0_RK
!    freqs_loop = [0:110]*0.1_RK
!    freqs_loop = [[-110:-90]*0.1_RK, [-8:8]*1.0_RK, [90:110]*0.1_RK]
!    freqs_loop = freqs_loop*1.0E-3_RK + freq01

!    phis_loop = [0.0_RK, pi] !SOF
    nPhis = 4 !Remember to clear files if you decrease this
    phis_loop = [0:nPhis-1]*2.0_RK*pi/nPhis !FOSOF
!    phis_loop = [0:nPhis-1]*2.0_RK*pi/nPhis*(1.0_RK-1.0E-3_RK) !Phases corresponding to offset freq error

    freqs_loop = [29609.458, 29610.583, 29611.708, 29612.833, 29613.958, 29614.458,&
                  29614.625, 29614.791, 29614.958, 29615.125, 29615.291, 29615.458,&
                  29615.625, 29615.791, 29615.958, 29616.125, 29616.291, 29616.458,&
                  29616.625, 29616.791, 29616.958, 29617.125, 29617.291, 29617.458,&
                  29617.625, 29617.791, 29617.958, 29618.125, 29618.291, 29618.458,&
                  29618.625, 29618.791, 29618.958, 29619.125, 29619.291, 29619.458,&
                  29619.958, 29621.083, 29622.208, 29623.333, 29624.458]*1.0E-3_RK
    nFreqs = SIZE(freqs_loop)
    f0Ind = NINT(nFreqs/2.0_RK)
    nPhis = SIZE(phis_loop)
    nPhis_overall = SIZE(phis_overall)
    nPhis_overall_float = REAL(nPhis_overall,RK)
    psi_vec = psi_vec/SQRT(SUM(ABS(psi_vec)**2)) !Normalise wavefunction

    allocate(B_background(nPhis_overall,3))
    B_background(1:nPhis_overall, 1) = 0.0_RK !Residual background B field (x,y,z)
    B_background(1:nPhis_overall, 2) = 0.0_RK
    B_background(1:nPhis_overall, 3) = 0.0_RK

	do parLoopInd=1,SIZE(loopArr) !Whole run loop for changing parameters
		loopVar = loopArr(parLoopInd)
!        Bz0 = loopVar
        powerLvl = loopVar

		allocate(psi_new(nPhis_overall,nStates), psi(nPhis_overall,nStates))
		allocate(B(nPhis_overall,3), B_Step(nPhis_overall,3), B_halfStep(nPhis_overall,3))

		do expNum=1,size(expTypes) !A/B experiment loop
			expType = expTypes(expNum)
			t_max = tT + tD + 20.0_RK
			centreT = t_max/2.0_RK

			do j=1,nPhis !Phis pulse loop
				write(loopStr, '(E25.17)') loopVar
!				firstSlash = 35!4
!				dividerLoc = index(subDir(firstSlash:), '\')+firstSlash
				resultsDir = resultsStem//'/Saturation effect at start of one pulse'!//trim(adjustl(subDir(dividerLoc:)))
				fullSaveDir = trim(adjustl(resultsDir))//'/'//trim(adjustl(loopStr))
				folStat = MAKEDIRQQ(trim(adjustl(resultsDir)))
				folStat = MAKEDIRQQ(trim(adjustl(fullSaveDir)))
				folStat = MAKEDIRQQ(trim(adjustl(fullSaveDir))//'/'//'results_'//expType)

				write(jstr,*) j
                open(20+j, file=trim(adjustl(fullSaveDir))//'/'//'results_'//expType//'/'//trim(adjustl(jstr))//'.csv', status='unknown')
				do i=1,nFreqs !Freqs loop
					freqRF = freqs_loop(i)
					tT = 300.0_RK !+ timeShiftScale*(freqRF-freq01)*1.0E3_RK
!					ampScale = 0.5_RK*SQRT(powerLvl/100.0_RK)*SQRT(1.0_RK+((freqRF-freq01)*10.0_RK*PSlope)) !*10 instead of /100 because GHz/MHz conversion
                    ampScale = 0.5_RK*SQRT(powerLvl/100.0_RK)!*SQRT(1.0_RK+((freqRF-freq01)*10.0_RK*PSlope)

					!Runge-Kutta initialising
					t = 0.0_RK
					do init_i=1,nPhis_overall
						psi_new(init_i, 1:nStates) = psi_vec
					end do
                    B_step(1:nPhis_overall, 1) = 0.0_RK
                    B_step(1:nPhis_overall, 2) = 0.0_RK
                    B_step(1:nPhis_overall, 3) = Bz0

					open(12, file=trim(adjustl(fullSaveDir))//'/'//'results_'//expType//'/'//'trace.csv', status='unknown')
!                            open(13, file=trim(adjustl(fullSaveDir))//'/'//'results_'//expType//'/'//'extras.csv', status='unknown')
					open(14, file=trim(adjustl(fullSaveDir))//'/'//'results_'//expType//'/'//'pops.csv', status='unknown')
					do while (t.LT.t_max) !Runge-Kutta loop

						if (t .GT. centreT) then
							if (expType .EQ. 'A') then
								phi = phis_loop(j)
!								phi_chirp = 0.0_RK
							else
								phi = 0.0_RK
!								phi_chirp = pi/2*((0.5+0.5*TANH((t-9)/chirpStretch)) - (0.5+0.5*TANH((t-11)/chirpStretch)) + &
!								 (0.5+0.5*TANH((t-(tD+9))/chirpStretch)) - (0.5+0.5*TANH((t-(tD+11))/chirpStretch)) + &
!                                 (0.5+0.5*TANH((t-(tT+9))/chirpStretch)) - (0.5+0.5*TANH((t-(tT+11))/chirpStretch)) + &
!                                 (0.5+0.5*TANH((t-(tT+tD+9))/chirpStretch)) - (0.5+0.5*TANH((t-(tT+tD+11))/chirpStretch)))
							end if
						else if (t.LT.centreT) then
							if (expType .EQ. 'B') then
								phi = phis_loop(j)
!								phi_chirp = 0.0_RK
							else
								phi = 0.0_RK
!								phi_chirp = pi/2*((0.5+0.5*TANH((t-9)/chirpStretch)) - (0.5+0.5*TANH((t-11)/chirpStretch)) + &
!								 (0.5+0.5*TANH((t-(tD+9))/chirpStretch)) - (0.5+0.5*TANH((t-(tD+11))/chirpStretch)) + &
!                                 (0.5+0.5*TANH((t-(tT+9))/chirpStretch)) - (0.5+0.5*TANH((t-(tT+11))/chirpStretch)) + &
!                                 (0.5+0.5*TANH((t-(tT+tD+9))/chirpStretch)) - (0.5+0.5*TANH((t-(tT+tD+11))/chirpStretch)))
							end if
						end if
						if (expType .EQ. 'A') then
                            phi_chirp = pi/2*(chirpScales(1)*((0.5+0.5*TANH((t-9)/chirpStretch)) - (0.5+0.5*TANH((t-11)/chirpStretch))) + &
                                     chirpScales(2)*((0.5+0.5*TANH((t-(tD+9))/chirpStretch)) - (0.5+0.5*TANH((t-(tD+11))/chirpStretch))) + &
                                     chirpScales(3)*((0.5+0.5*TANH((t-(tT+9))/chirpStretch)) - (0.5+0.5*TANH((t-(tT+11))/chirpStretch))) + &
                                     chirpScales(4)*((0.5+0.5*TANH((t-(tT+tD+9))/chirpStretch)) - (0.5+0.5*TANH((t-(tT+tD+11))/chirpStretch))))
                        else
                            phi_chirp = pi/2*(chirpScales(3)*((0.5+0.5*TANH((t-9)/chirpStretch)) - (0.5+0.5*TANH((t-11)/chirpStretch))) + &
                                     chirpScales(4)*((0.5+0.5*TANH((t-(tD+9))/chirpStretch)) - (0.5+0.5*TANH((t-(tD+11))/chirpStretch))) + &
                                     chirpScales(1)*((0.5+0.5*TANH((t-(tT+9))/chirpStretch)) - (0.5+0.5*TANH((t-(tT+11))/chirpStretch))) + &
                                     chirpScales(2)*((0.5+0.5*TANH((t-(tT+tD+9))/chirpStretch)) - (0.5+0.5*TANH((t-(tT+tD+11))/chirpStretch))))
                        end if
!                        phi_chirp = ((-0.4_RK*(t-8.5_RK) + 1.5_RK)*(t.GT.0.0_RK)*(1+tanh((-t+12.0_RK)*8.0_RK))+ &
!                                     (-0.4_RK*(t-308.5_RK) + 1.5_RK)*(t.GT.300.0_RK)*(1+tanh((-t+312.0_RK)*8.0_RK)))*0.5_RK
                        phi_chirp = ((-0.4_RK*(t-5.5_RK) + 1.5_RK)*(t.GT.0.0_RK)*(1+TANH((-t+9.0_RK)*8.0_RK))+ &
                                     (-0.4_RK*(t-305.5_RK) + 1.5_RK)*(t.GT.300.0_RK)*(1+TANH((-t+309.0_RK)*8.0_RK)))*0.5_RK
						phis = phis_overall + phi + phi_chirp*chirpScale

						!z-polarised microwaves with some y component, and z-aligned DC B field. Use previous stepped value to save on calculations
						psi = psi_new
						B = B_step
						B_halfStep(1:nPhis_overall, 1) = 0.0_RK
						B_halfStep(1:nPhis_overall, 2) = getB_mw(t+h/2.0_RK)*SIND(theta_MW)
						B_halfStep(1:nPhis_overall, 3) = Bz0 + getB_mw(t+h/2.0_RK)*COSD(theta_MW)
						B_step(1:nPhis_overall, 1) = 0.0_RK
						B_step(1:nPhis_overall, 2) = getB_mw(t+h)*SIND(theta_MW)
						B_step(1:nPhis_overall, 3) = Bz0 + getB_mw(t+h)*COSD(theta_MW)

						k1 = schrodStep(EJ0, EJ1, EJ2, psi, B+B_background)
						k2 = schrodStep(EJ0, EJ1, EJ2, psi+h*k1/2.0_RK, B_halfStep+B_background)
						k3 = schrodStep(EJ0, EJ1, EJ2, psi+h*k2/2.0_RK, B_halfStep+B_background)
						k4 = schrodStep(EJ0, EJ1, EJ2, psi+h*k3, B_step+B_background)
						k = (k1 + 2.0_RK*k2 + 2.0_RK*k3 + k4)/6.0_RK

						psi_new = psi + h*k
						t = t+h

						if((j.EQ.nPhis).AND.(i.EQ.f0Ind).AND.(mod(t,saveStep).LT.h)) then
							write(12,WriteFmtStr2) t, ",", amp, ",", phis(2)
							write(14,WriteFmtStr3) t, ",", sum(abs(psi_new(1:nPhis_overall,1))**2)/nPhis_overall_float, &
							",", sum(abs(psi_new(1:nPhis_overall,2))**2)/nPhis_overall_float, &
							",", sum(abs(psi_new(1:nPhis_overall,3))**2)/nPhis_overall_float, &
							",", sum(abs(psi_new(1:nPhis_overall,4))**2)/nPhis_overall_float
						end if

					end do !Runge-Kutta time loop
					close(12)
					close(14)
					write(20+j,WriteFmtStr) freqRF, ",", sum(abs(psi_new(1:nPhis_overall,1))**2)/nPhis_overall_float
					write(loopSizeStr, '(I2,A,I2)') parLoopInd,'/',size(loopArr)
					write(*,*) trim(loopSizeStr)//' ', expType, j, freqRF
				end do !Freqs loop
                close(20+j)
			end do !Phis pulse loop
			deallocate(phis)
		end do !A/B experiment loop
		deallocate(B, B_halfStep, B_Step, psi, psi_new)
	end do !Whole run loop for changing parameters

contains
    function getB_mw(t)
    !Gives two-pulse profile, scaled
        real(RK), intent(in) :: t
        real(RK), dimension(nPhis_overall) :: getB_mw
        real(RK) :: t_shift

        t_shift = timeShiftScale*(freqRF-freq01)*1.0E3_RK

!        if (t .LT. 10.0_RK) then
!			amp = 0.0_RK
!		else if (t .LT. 10.0_RK+tD) then
!!            amp = 1.0_RK
!            if (expType .EQ. 'A') then
!                amp = SQRT(1.0_RK-PwrImbalance/200.0_RK)*SQRT(1.0_RK+((freqRF-freq01)*10.0_RK*PSlope))
!            else if (expType .EQ. 'B') then
!                amp = SQRT(1.0_RK+PwrImbalance/200.0_RK)
!            end if
!		else if (t .LT. 10.0_RK+tT) then
!			amp = 0.0_RK
!		else if (t .LT. 10.0_RK+tT+tD) then
!!			amp = 1.0_RK
!            if (expType .EQ. 'A') then
!                amp = SQRT(1.0_RK+PwrImbalance/200.0_RK)
!            else if (expType .EQ. 'B') then
!                amp = SQRT(1.0_RK-PwrImbalance/200.0_RK)*SQRT(1.0_RK+((freqRF-freq01)*10.0_RK*PSlope))
!            end if
!		else if (t .LE. t_max) then
!			amp = 0.0_RK
!        end if

!        if (expType .EQ. 'A') then
!            if (t .LT. 10.0_RK+t_shift) then
!                amp = 0.0_RK
!            else if (t .LT. 10.0_RK+t_shift+tD) then
!                amp = 1.0_RK
!            else if (t .LT. 10.0_RK+tT) then
!                amp = 0.0_RK
!            else if (t .LT. 10.0_RK+tT+tD) then
!                amp = 1.0_RK
!            else
!                amp = 0.0_RK
!            end if
!        else if (expType .EQ. 'B') then
!            if (t .LT. 10.0_RK) then
!                amp = 0.0_RK
!            else if (t .LT. 10.0_RK+tD) then
!                amp = 1.0_RK
!            else if (t .LT. 10.0_RK+tT+t_shift) then
!                amp = 0.0_RK
!            else if (t .LT. 10.0_RK+tT+tD+t_shift) then
!                amp = 1.0_RK
!            else
!                amp = 0.0_RK
!            end if
!        end if
        if (expType .EQ. 'A') then
            amp = (TANH(t-10.0_RK)-TANH(t-tD-10.0_RK)+TANH(t-tT-10.0_RK)-TANH(t-tT-tD-10.0_RK))/2 + &
                   0.008_RK*powerLvl*(EXP(-(t-9.0_RK)**2.0_RK/2.0_RK))!+EXP(-(t-309.0_RK)**2.0_RK/2.0_RK))
        else if (expType .EQ. 'B') then
            amp = (TANH(t-10.0_RK)-TANH(t-tD-10.0_RK)+TANH(t-tT-10.0_RK)-TANH(t-tT-tD-10.0_RK))/2 + &
                   0.008_RK*powerLvl*EXP(-(t-309.0_RK)**2.0_RK/2.0_RK)
        end if
        amp = amp*ampScale
        getB_mw = amp*SIN(2*pi*freqRF*t+phis)! + ampScale*SQRT(1.0E-4_RK)*(SIN(2*pi*freqRF*t+phis_overall))!+SIN(2*pi*freqRF*t+phis_loop(j))) !2nd term simulates bad switches
    end function

    function schrodStep(EJ0, EJ1, EJ2, psi, B)
    !Calculate the 9-level Schrodinger eqn
        real(RK), intent(in) :: EJ0, EJ1, EJ2
        real(RK), dimension(nPhis_overall, 3), intent(in) :: B
		real(RK), dimension(nPhis_overall) :: Bx, By, Bz
        complex(RK), dimension(nPhis_overall, nStates), intent(in) :: psi
        complex(RK), dimension(nPhis_overall, nStates) :: schrodStep
        complex(RK), dimension(nPhis_overall) :: psiJ0m0, psiJ1mn1, psiJ1m0, psiJ1m1, psiJ2mn2, psiJ2mn1, psiJ2m0, psiJ2m1, psiJ2m2

		Bx = B(1:nPhis_overall, 1)
		By = B(1:nPhis_overall, 2)
		Bz = B(1:nPhis_overall, 3)

		psiJ0m0 = psi(1:nPhis_overall, 1)
		psiJ1mn1 = psi(1:nPhis_overall, 2)
		psiJ1m0 = psi(1:nPhis_overall, 3)
		psiJ1m1 = psi(1:nPhis_overall, 4)
		psiJ2mn2 = psi(1:nPhis_overall, 5)
		psiJ2mn1 = psi(1:nPhis_overall, 6)
		psiJ2m0 = psi(1:nPhis_overall, 7)
		psiJ2m1 = psi(1:nPhis_overall, 8)
		psiJ2m2 = psi(1:nPhis_overall, 9)

        schrodStep(1:nPhis_overall, 1) = (xi*psiJ1mn1*(Bx-ci*By)/SQRT(3.0_RK) - xi*psiJ1m1*(Bx+ci*By)/SQRT(3.0_RK) &
        + SQRT(2.0_RK/3.0_RK)*xi*psiJ1m0*Bz + psiJ0m0*EJ0)*-ci/hbar

        schrodStep(1:nPhis_overall, 2) = (0.5_RK*xi*psiJ2mn2*(Bx-ci*By) + xi*psiJ0m0*(Bx+ci*By)/SQRT(3.0_RK) &
        + chi*psiJ1m0*(Bx+ci*By)/SQRT(8.0_RK) - xi*psiJ2m0*(Bx+ci*By)/SQRT(24.0_RK) + 0.5_RK*xi*psiJ2mn1*Bz &
        + psiJ1mn1*(EJ1-chi*Bz/2.0_RK))*-ci/hbar

		schrodStep(1:nPhis_overall, 3) = (chi*psiJ1mn1*(Bx-ci*By)/SQRT(8.0_RK) + xi*psiJ2mn1*(Bx-ci*By)/SQRT(8.0_RK) &
		+ chi*psiJ1m1*(Bx+ci*By)/SQRT(8.0_RK) - xi*psiJ2m1*(Bx+ci*By)/SQRT(8.0_RK) + SQRT(2.0_RK/3.0_RK)*xi*psiJ0m0*Bz &
		+ xi*psiJ2m0*Bz/SQRT(3.0_RK) + psiJ1m0*EJ1)*-ci/hbar

		schrodStep(1:nPhis_overall, 4) = (-xi*psiJ0m0*(Bx-ci*By)/SQRT(3.0_RK) + chi*psiJ1m0*(Bx-ci*By)/SQRT(8.0_RK) &
		+ xi*psiJ2m0*(Bx-ci*By)/SQRT(24.0_RK) - 0.5_RK*xi*psiJ2m2*(Bx+ci*By) + 0.5_RK*xi*psiJ2m1*Bz &
		+ psiJ1m1*(EJ1+chi*Bz/2.0_RK))*-ci/hbar

		schrodStep(1:nPhis_overall, 5) = (0.5_RK*xi*psiJ1mn1*(Bx+ci*By) + 0.5_RK*chi*psiJ2mn1*(Bx+ci*By) &
		+ psiJ2mn2*(EJ2-chi*Bz))*-ci/hbar

		schrodStep(1:nPhis_overall, 6) = (0.5_RK*chi*psiJ2mn2*(Bx-ci*By) + xi*psiJ1m0*(Bx+ci*By)/SQRT(8.0_RK) &
		+ SQRT(3.0_RK/8.0_RK)*chi*psiJ2m0*(Bx+ci*By) + 0.5*xi*psiJ1mn1*Bz + psiJ2mn1*(EJ2-chi*Bz/2.0_RK))*-ci/hbar

		schrodStep(1:nPhis_overall, 7) = (-xi*psiJ1mn1*(Bx-ci*By)/SQRT(24.0_RK) &
		+ SQRT(3.0_RK/8.0_RK)*chi*psiJ2mn1*(Bx-ci*By) + xi*psiJ1m1*(Bx+ci*By)/SQRT(24.0_RK) &
		+ SQRT(3.0_RK/8.0_RK)*chi*psiJ2m1*(Bx+ci*By) + xi*psiJ1m0*Bz/SQRT(3.0_RK) + psiJ2m0*EJ2)*-ci/hbar

		schrodStep(1:nPhis_overall, 8) = (-xi*psiJ1m0*(Bx-ci*By)/SQRT(8.0_RK) + SQRT(3.0_RK/8.0_RK)*chi*psiJ2m0*(Bx-ci*By) &
		+ 0.5_RK*chi*psiJ2m2*(Bx+ci*By) + 0.5_RK*xi*psiJ1m1*Bz + psiJ2m1*(EJ2 + chi*Bz/2.0_RK))*-ci/hbar

		schrodStep(1:nPhis_overall, 9) = (-0.5_RK*xi*psiJ1m1*(Bx-ci*By) + 0.5_RK*chi*psiJ2m1*(Bx-ci*By) &
        + psiJ2m2*(EJ2 + chi*Bz))*-ci/hbar

    end function


end program

