program main
		
		parameter (NJ=501,NI=101,NK=4,NL=2,N=10000000)
		double precision EP, RV(NI), ALPHA, BETA, curr(NI)
		double precision Q, M, T, KB, SUM(NI), tol(3), pi, norm, vol
		double precision vx, vy, vz, vx0, vy0, vz0, vth, chi(NI,NJ,NK,NL)
		double precision vxl, vxu, vyl, vyu, vzl, vzu, vr, del(3), chi0
		double precision probx, proby, probz, prob, kappa, ncur, psi
		REAL ETIME, ELAPSED(2)	
		double precision dx, dy, dz, vcross
		
C-------------------------------------
C  Initialize variables
C-------------------------------------

		  ISEED = -1237
		  iseeda = -1489
		  
		  IN = 10

	      Q = 1.60217646E-19
	      M = 15.9994*1.67262158E-27
		  KB = 1.3806503E-23
		  
C-------------------------------------
C  Load transmission coefficients
C-------------------------------------
	  
  		  open(2,file="debug-sq.dat")
          idum=0
          nn = ni*nj*nk*nl
          do idum = 1, nn
             READ(2,*) i,j,k,l, chi0
             chi(i,j,k,l) = chi0
          enddo
      
          close(2)
          
		  open(1,file='iv-randomgen3-key.dat')
          open(3,file="iv-randomgen3.dat")
          
          write(1,*) IN



C-------------------------------------
C  Begin Big Loop
C-------------------------------------
	  
		  do jdum = 1, IN

			 write(6,102) jdum

		     vx0 = 1600.*(ran1(iseeda)-0.5)
		     vy0 = 1600.*(ran1(iseeda)-0.5)
		     vz0 =  7000. + 2000.*ran1(iseeda)
		     
		     vcross = sqrt(vx0**2 + vy0**2)
		     
			 T = 1000. + 4000.*ran1(iseeda)
			 
			 write(1,102) jdum
			 write(1,*) vx0, vy0, vz0
			 write(1,*) vz0, vcross, T
			 
		     vth = sqrt(2.*kb*T/m)
		  	  
		     tol(1) = 1.e-3
		     tol(2) = 1.e-3
		     tol(3) = 1.e-3
	  
		     PI = ACOS(-1.)
		   
		     norm = 0.
		  
		     do i = 1, ni
		        RV(i) = 0.1*(i-1)
		        sum(i) = 0.
		     enddo	  

C------------------------------------------------
C  Load center values
C------------------------------------------------
	
		     do i = 1, 3
		        del(i) = vth*sqrt(-log(tol(i)))
		     enddo
		 
		     vxl = vx0 - del(1)
		     vxu = vx0 + del(1)
		     dx = vxu-vxl
		
		     vyl = vy0 - del(2)
		     vyu = vy0 + del(2)
		     dy = vyu-vyl
		
		     vzl = 0.
		     vzu = vz0 + del(3)
		     dz = vzu-vzl
		     
		  ISEED = -1237 - 10*jdum
		     
	
C------------------------------------------------
C  Begin Little Loop
C------------------------------------------------

		     do idum = 1, N
	 
C------------------------------------------------
C  Generate random coordinates and integrate
C------------------------------------------------

	   	        vx = vxl + dx*ran2(iseed)
		        vy = vyl + dy*ran2(iseed)
		        vz = vzl + dz*ran2(iseed)
		   
		        EP = 0.5*m*vz**2/q
		        alpha = (180./pi) * atan(sqrt(vx**2+vy**2)/vz)
		        beta = (180./pi) * atan(vx/vy)
		        if (beta.lt.0.) beta = beta + 90.
		   
		        if(alpha.lt.45.) then
		           j = int(50*EP) + 1
		           if (j.gt.nj) j = nj
		           k = int(alpha/5.-0.5) + 1
		           if (k.gt.nk) k = nk
		           if (k.eq.0) k=1
		           l = int(beta/45.) + 1
		           if (l.gt.nl) l = l-2
		   
		           probx = exp(-((vx-vx0)/vth)**2)
		           proby = exp(-((vy-vy0)/vth)**2)
		           probz = exp(-((vz-vz0)/vth)**2)
		   
		           prob = probx*proby*probz
		   
		           norm = norm + prob
		      
		           psi = vz*prob

		           do i = 1, ni
		              sum(i) = sum(i) + chi(i,j,k,l)*psi
		           enddo
		      
		        endif
		     enddo
		
		     vol = norm*dx*dy*dz/n
		     vol = vol /(sqrt(pi)*vth)**3.

C------------------------------------------------
C  Normalize Currents
C------------------------------------------------
		
		     do i = 1, ni
		        curr(i) = sum(i)/sum(1)
		     enddo
		
C------------------------------------------------
C  Output to file
C------------------------------------------------
		
             write(3,*) 'VARIABLES = "Voltage","Whipple","MC"'
             write(3,101) jdum, vol
             do i = 1, ni
                vr = sqrt(2.*q*RV(i)/m)
	            KAPPA = (VZ0 - vr)/vth
                ncur = ERF(KAPPA)+EXP(-KAPPA**2.)/(sqrt(PI)*VZ0/vth)
                ncur = 0.5*(1.+ncur)
             
                write(3,100) RV(i), ncur, curr(i)
             enddo
          enddo
          close(1)
          close(3)
            
 100    format(f4.1,3x,f6.4,3x,f6.4)	
 101    format('ZONE T="Case ',I4,',Vol=',F6.4,'", I=101')
 102    format('Case ',I4)

		
        if(etime(elapsed).lt.60.) then
           WRITE(6,110) ETIME(ELAPSED)
        elseif(etime(elapsed).lt.3600.) then
           WRITE(6,111) ETIME(ELAPSED)/60.
        elseif(etime(elapsed).lt.86400.) then
           WRITE(6,112) ETIME(ELAPSED)/3600.
		else
           WRITE(6,113) ETIME(ELAPSED)/86400.
        endif

 110    format(' Elapsed Time is ', f6.2,' seconds')        
 111    format(' Elapsed Time is ', f6.2,' minutes')        
 112    format(' Elapsed Time is ', f6.2,' hours')        
 113    format(' Elapsed Time is ', f6.2,' days')        
		
		
		end
************************************************************

		include "ran2.f"
		include "ran1.f"