module initialconfignewmodule
contains
subroutine initialconfignew(N,length,seed,rx,ry,rz,vx,vy,vz)

    !use velrescalemodule
	implicit none
	integer, intent(IN) :: N,seed 	!global parameters
	double precision, intent(IN) :: length
	double precision, dimension(:), intent(OUT) :: rx,ry,rz,vx,vy,vz   !global parameters
	double precision q,halfbox,v,r,s                  !,Vmax=0.2	!local parameters
	integer :: i,j,k,g,p,b,m,ab 
    double precision, dimension(4) :: x1,y1,z1

    halfbox=length/2.0             !for shifiting axis to (h,h,h), @ centre of cube
	p=(N**(1./3.))		!fixing number of particles to place in each direction
	!q=L*1.0/p			!cal. equal spacing for p particles in L space
	g=1


!!assigning initial coordinates to particles

!!! electron coordinates
s=2.5
r=3.0
v=4.5
    do i=1,92

        rx(i)=21.5*length*RAND()
        ry(i)=21.5*length*RAND()
        rz(i)=21.5*length*RAND()

        s=s+0.4
        r=r+0.3
        v=v-0.4

    end do


!!! proton and neutron coordinates
s=2.0
r=1.5
v=3.0
    do i=93,330

        rx(i)=5.0+3.5*RAND()      !placing at centre of box
        ry(i)=5.0+3.5*RAND()
        rz(i)=5.0+3.2*RAND()

        s=s+0.04
        r=r+0.02
        v=v-0.02

    end do



	!giving initial velocities acc. to boltzmann normal distribution with scaling temp T
	b=seed
	do i=1,N
        b=seed
		call boltzmann(i,b,v)   !i,b serves as a seed
        vx(i)=v
        b=b+8223
		call boltzmann(i,b,v)
        vy(i)=v
        b=b+9334321
        call boltzmann(i,b,v)
        vz(i)=v

	end do


!!! velocity of proton n neutron to zero
    do i=9,24
        vx(i)=0.0
        vy(i)=0.0
        vz(i)=0.0
    end do

!!! end of velocity to zero



 !call velrescale(N,vx,vy,vz,Temp)



end subroutine initialconfignew




    subroutine boltzmann(i,b,v)
        implicit none
        integer, intent(IN) :: i
        integer, intent(inout) :: b
        double precision, intent(out) :: v
        double precision a,x,y
        integer u,w

        u=i+b
        w=i+b
		do

			x=rand((i+207*i)*u)
			y=rand((i+223*i)*w)
			x=2.0*x-1.0                 !interval of -1 to 1
			y=2.0*y-1.0                 !interval of -1 to 1
			a=(x*x)+(y*y)       !a should be in between 0 to 1
			if(a>0.0.AND.a<1.0)exit
            u=int(x*10000)
            w=int(y*10000)

		end do
		v=x*sqrt(-2.0*log(a)/a)
	end subroutine boltzmann
		

end module initialconfignewmodule
	  
