MODULE diffusion
	USE interpol
	USE spalation
	USE gauss
	USE arrays
	IMPLICIT NONE
	CONTAINS
	SUBROUTINE diffusions(xsec,disc_r,disc_v,cooling_E,cooling_b,nuc_density,events,density)
		IMPLICIT NONE
		REAL(8),INTENT(in)::xsec(min_Z:max_Z,min_A:max_A,min_cell_E:max_cell_E,min_Z:max_Z,min_A:max_A,min_cell_E:max_cell_E)
		REAL(8),INTENT(in)::disc_r(1:n),disc_v(1:n),cooling_E(1:m),cooling_b(1:m)
		INTEGER,INTENT(in out)::events(min_Z:max_Z,min_A:max_A,min_cell_E:max_cell_E)
		REAL(8),INTENT(in out)::nuc_density(min_Z:max_Z,min_A:max_A,min_cell_E:max_cell_E)
		REAL(4),INTENT(in out)::density(-cell_XY:cell_XY,-cell_XY:cell_XY,-cell_Z:cell_Z,min_cell_E:max_cell_E)
		INTEGER::j,e,ii
		REAL(8)::initial_energy,v,coords(3),omega,weight
		weight=1							!****************check this**********
		DO j=1,samples
			CALL source(coords,omega)

			SELECT CASE(calculate)
				CASE('S/P')
					CALL initial_particle()
			END SELECT 
			CALL initial(e,initial_energy,v)
					
			CALL particle_diffusion(initial_energy,omega,xsec,disc_r,disc_v,cooling_E,cooling_b,&
									&e,nuc_density,events,density,v,coords,weight)
		END DO
	END SUBROUTINE diffusions

	SUBROUTINE particle_diffusion(initial_energy,omega,xsec,disc_r,disc_v,cooling_E,cooling_b,&
								&e,nuc_density,events,density,v,coords,weight)
		IMPLICIT NONE
		REAL(8),INTENT(in)::initial_energy,omega
		REAL(8),INTENT(in)::xsec(min_Z:max_Z,min_A:max_A,min_cell_E:max_cell_E,min_Z:max_Z,min_A:max_A,min_cell_E:max_cell_E)
		REAL(8),INTENT(in)::disc_r(1:n),disc_v(1:n),cooling_E(1:m),cooling_b(1:m)
		INTEGER,INTENT(in out)::e,events(min_Z:max_Z,min_A:max_A,min_cell_E:max_cell_E)
		REAL(8),INTENT(in out)::nuc_density(min_Z:max_Z,min_A:max_A,min_cell_E:max_cell_E)
		REAL(4),INTENT(in out)::density(-cell_XY:cell_XY,-cell_XY:cell_XY,-cell_Z:cell_Z,min_cell_E:max_cell_E)
		REAL(8),INTENT(in out)::v,coords(3),weight
		INTEGER::e_sec
		REAL(8)::energy,dt,l,g,pri_energy,pri_coords(3),q
		LOGICAL::break
		energy=initial_energy
		g=0
		break=.false.
		DO
			IF(SQRT(coords(1)**2+coords(2)**2)>r_galaxy.OR.ABS(coords(3))>z_galaxy)EXIT
			IF(e<min_cell_E.OR.e>max_cell_E)EXIT
			SELECT CASE(calculate)
				CASE('S/P')
					SELECT CASE(model)
						CASE('spiral arms')
							IF(SQRT((ABS(coords(1))-8.5)**2+coords(2)**2)>z_galaxy*5)EXIT
					END SELECT 
			END SELECT
			
			SELECT CASE(calculate)
				CASE('S/P')
					CALL particle_spalation(e,initial_energy,energy,v,g,coords,weight,xsec,nuc_density,events,e_sec,dt,l,break)
				
				CASE('galaxy')
					CALL galaxy_density(e,initial_energy,v,coords,energy,density,dt,l)
			END SELECT 

			CALL grammage(g,dt,v,coords)

			CALL energy_losses(g,dt,cooling_E,cooling_b,energy,v,e)

			CALL step(dt,l,omega,disc_r,disc_v,coords)
			
			SELECT CASE(plot)
				CASE('B/C','Iron')
					IF(break)EXIT
				CASE('PAMELA')
					IF(break)THEN
						e=e_sec
						pri_energy=energy
						pri_coords=coords
						CALL secondary_particle_diffusion(initial_energy,omega,xsec,disc_r,disc_v,cooling_E,cooling_b,&
														&e,nuc_density,events,density,v,coords,weight)
						Z=1
						A=1
						break=.false.
						energy=pri_energy-10**(e_sec*energy_grid)			!***************check this***************
						q=log10(energy)
						IF(ABS((q/energy_grid)-int(q/energy_grid))<0.5)THEN
							e=int(q/energy_grid)
						ELSE
							IF(q>0)THEN
								e=int(q/energy_grid)+1
							ELSE
								e=int(q/energy_grid)-1
							END IF
						END IF
						v=SQRT(1-1/(1+energy)**2)
						coords=pri_coords
					END IF
				END SELECT
		END DO
	END SUBROUTINE particle_diffusion
	
	SUBROUTINE secondary_particle_diffusion(initial_energy,omega,xsec,disc_r,disc_v,cooling_E,cooling_b,&
								&e,nuc_density,events,density,v,coords,weight)
		IMPLICIT NONE
		REAL(8),INTENT(in)::initial_energy,omega
		REAL(8),INTENT(in)::xsec(min_Z:max_Z,min_A:max_A,min_cell_E:max_cell_E,min_Z:max_Z,min_A:max_A,min_cell_E:max_cell_E)
		REAL(8),INTENT(in)::disc_r(1:n),disc_v(1:n),cooling_E(1:m),cooling_b(1:m)
		INTEGER,INTENT(in out)::e,events(min_Z:max_Z,min_A:max_A,min_cell_E:max_cell_E)
		REAL(8),INTENT(in out)::nuc_density(min_Z:max_Z,min_A:max_A,min_cell_E:max_cell_E)
		REAL(4),INTENT(in out)::density(-cell_XY:cell_XY,-cell_XY:cell_XY,-cell_Z:cell_Z,min_cell_E:max_cell_E)
		REAL(8),INTENT(in out)::v,coords(3),weight
		INTEGER::e_sec
		REAL(8)::energy,dt,l,g
		LOGICAL::break
		energy=10**(e*energy_grid)
		v=1
		g=0
		break=.false.
		DO
			IF(SQRT(coords(1)**2+coords(2)**2)>r_galaxy.OR.ABS(coords(3))>z_galaxy)EXIT
			IF(e<min_cell_E.OR.e>max_cell_E)EXIT
			SELECT CASE(calculate)
				CASE('S/P')
					SELECT CASE(model)
						CASE('spiral arms')
							IF(SQRT((ABS(coords(1))-8.5)**2+coords(2)**2)>z_galaxy*5)EXIT
					END SELECT 
			END SELECT
			
			SELECT CASE(calculate)
				CASE('S/P')
					CALL particle_spalation(e,initial_energy,energy,v,g,coords,weight,xsec,nuc_density,events,e_sec,dt,l,break)
				
				CASE('galaxy')
					CALL galaxy_density(e,initial_energy,v,coords,energy,density,dt,l)
			END SELECT 

			CALL grammage(g,dt,v,coords)

			CALL energy_losses(g,dt,cooling_E,cooling_b,energy,v,e)

			CALL step(dt,l,omega,disc_r,disc_v,coords)
		END DO
	END SUBROUTINE secondary_particle_diffusion
	
	SUBROUTINE galaxy_density(e,initial_energy,v,coords,energy,density,dt,l)
		IMPLICIT NONE
		INTEGER,INTENT(in)::e
		REAL(8),INTENT(in)::initial_energy,v,coords(3),energy
		REAL(4),INTENT(in out)::density(-cell_XY:cell_XY,-cell_XY:cell_XY,-cell_Z:cell_Z,min_cell_E:max_cell_E)
		REAL(8),INTENT(in out)::dt,l
		INTEGER::j,cell(3)
		REAL(8)::p,D,t
		DO j=1,2
			IF((ABS(coords(j)/grid-int(coords(j)/grid)))<0.5)THEN
				cell(j)=coords(j)/grid
			ELSE
				IF(coords(j)<0)THEN
					cell(j)=coords(j)/grid-1
				ELSE
					cell(j)=coords(j)/grid+1
				END IF
			END IF
		END DO
		IF((ABS(coords(3)/(grid/10)-int(coords(3)/(grid/10))))<0.5)THEN
			cell(3)=coords(3)/(grid/10)
		ELSE
			IF(coords(3)<0)THEN
				cell(3)=coords(3)/(grid/10)-1
			ELSE
				cell(3)=coords(3)/(grid/10)+1
			END IF
		END IF
		SELECT CASE(plot)
			CASE('nuclei')
				Z=1
				A=2
				p=A*SQRT(energy**2+2*energy)
				D=D0*v*(p/(3*Z))**(delta)
			CASE('protons')
				Z=1
				A=1
				p=A*SQRT(energy**2+2*energy)
				D=D0*v*(p/(3*Z))**(delta)
			CASE('electrons')				!******add PAMELA******
				D=D0*v*energy**(delta)
		END SELECT 
		t=z_galaxy**2/(2*D)
		dt=0.01*t
		l=SQRT(2*D*dt)
		density(cell(1),cell(2),cell(3),e)=density(cell(1),cell(2),cell(3),e)+dt*(initial_energy**(-(1.7-delta)))
	END SUBROUTINE galaxy_density

	SUBROUTINE grammage(g,dt,v,coords)
		IMPLICIT NONE
		REAL(8),INTENT(in)::dt,v,coords(3)
		REAL(8),INTENT(out)::g
		REAL(8)::r,height,n_ISM(3)
		
			r=SQRT(coords(1)**2+coords(2)**2)
			
			height=ABS(coords(3))
			
			n_ISM(1)=0.015*EXP(-height/0.07)+0.025*EXP(-height/0.9)					!see Ferriere 2001 for ISM distribution data
			
			n_ISM(2)=0.57*(0.7*EXP(-(height/0.127)**2)+0.19*EXP(-(height/0.318)**2)+&
				&0.11*EXP(-height/0.403))*EXP(-MAX(0.,r-13.5)/4.)*MIN(1.,r/3.5)
				
			n_ISM(3)=0.53*EXP(-(height/0.071)**2)*EXP((8-MAX(4.,r))/4.)*MIN(1.,r/3.5)
			
			g=1.4*dt*megayear_sec*v*v_light_v_cm*(n_ISM(1)+n_ISM(2)+n_ISM(3))/(cm2_mb)
			
			IF(SQRT((coords(1)-8.5)**2+coords(2)**2+coords(3)**2)<0.1.OR.&
				&SQRT((coords(1)-8.35)**2+coords(2)**2+coords(3)**2)<0.1.OR.&
				&SQRT((coords(1)-8.2)**2+coords(2)**2+coords(3)**2)<0.15)g=0.1*g
	END SUBROUTINE grammage

	SUBROUTINE energy_losses(g,dt,cooling_E,cooling_b,energy,v,e)
		IMPLICIT NONE
		REAL(8),INTENT(in)::g,dt,cooling_E(1:m),cooling_b(1:m)
		REAL(8),INTENT(in out)::energy,v
		INTEGER,INTENT(out)::e
		REAL(8)::b,q
				SELECT CASE(plot)
					CASE('B/C','Iron','nuclei','protons')
						energy=energy-(g*cm2_mb*erg_GeV*4*pi*(r_e**2)*m_e*((v_light_v_cm/v)**2)*((Z**2)/A)*45/A)
						v=SQRT(1-1/(1+energy)**2)
						SELECT CASE(plot)
							CASE('B/C')
								half_t_Be10_rel=half_t_Be10/SQRT(1-v**2)
								half_t_Al26_rel=half_t_Al26/SQRT(1-v**2)
							CASE('Iron')
								half_t_Ni59_rel=half_t_Ni59/SQRT(1-v**2)
								half_t_Fe60_rel=half_t_Fe60/SQRT(1-v**2)
								half_t_Mn53_rel=half_t_Mn53/SQRT(1-v**2)
						END SELECT 
					CASE('PAMELA')
						IF(A==1)THEN
							energy=energy-(1.3*g*cm2_mb*erg_GeV*4*pi*(r_e**2)*m_e*((v_light_v_cm/v)**2)*((Z**2)/A)*45/A)
						ELSE IF(A==0)THEN
							b=interp(m,cooling_E,cooling_b,energy)
							energy=energy-(dt*megayear_sec*b*(10**(-16.))*(energy**2))
						END IF
					CASE('electrons')
						b=interp(m,cooling_E,cooling_b,energy)
						energy=energy-(dt*megayear_sec*b*(10**(-16.))*(energy**2))
				END SELECT 
				q=log10(energy)
				IF(ABS((q/energy_grid)-int(q/energy_grid))<0.5)THEN
					e=int(q/energy_grid)
				ELSE
					IF(q>0)THEN
						e=int(q/energy_grid)+1
					ELSE
						e=int(q/energy_grid)-1
					END IF
				END IF
	END SUBROUTINE energy_losses

	SUBROUTINE initial(e,initial_energy,v)
		IMPLICIT NONE
		INTEGER,INTENT(out)::e
		REAL(8),INTENT(out)::initial_energy,v
		REAL(8)::q
		CALL random_number(q)
		q=q*(max_cell_E-min_cell_E)*energy_grid+min_cell_E*energy_grid
		IF(ABS(((q/energy_grid)-int(q/energy_grid)/1))<0.5)THEN
			e=int(q/energy_grid)
		ELSE
			IF(q>0)THEN
				e=int(q/energy_grid)+1
			ELSE
				e=int(q/energy_grid)-1
			END IF
		END IF
		initial_energy=10**(q)
		IF(initial_energy<=500)THEN
			beta=1.75
			alpha=1
		ELSE
			beta=1.65
			alpha=500**(1.65-1.75)
		END IF
		SELECT CASE(plot)
			CASE('B/C','Iron','nuclei','protons')
				v=SQRT(1-1/(1+initial_energy)**2)
			CASE('PAMELA')
				IF(A==1)THEN
					v=SQRT(1-1/(1+initial_energy)**2)
				ELSE IF(A==0)THEN
					v=1
				END IF
			CASE('electrons')
				v=1
		END SELECT 
	END SUBROUTINE initial

	SUBROUTINE initial_particle()
		IMPLICIT NONE
		REAL(8)::q
		CALL random_number(q)
		SELECT CASE(plot)
			CASE('B/C')
				IF(q<per_Si*0.031)THEN
					Z=14
					A=30
				ELSE IF(q<per_Si*0.0777)THEN
					Z=14
					A=29
				ELSE IF(q<per_Si)THEN
					Z=14
					A=28
				ELSE IF(q<per_Si+per_Mg*0.1101)THEN
					Z=12
					A=26
				ELSE IF(q<per_Si+per_Mg*0.2101)THEN
					Z=12
					A=25
				ELSE IF(q<per_Si+per_Mg)THEN
					Z=12
					A=24
				ELSE IF(q<per_Si+per_Mg+per_Ne*0.0925)THEN
					Z=10
					A=22
				ELSE IF(q<per_Si+per_Mg+per_Ne*0.0952)THEN
					Z=10
					A=21
				ELSE IF(q<per_Si+per_Mg+per_Ne)THEN
					Z=10
					A=20
				ELSE IF(q<per_Si+per_Mg+per_Ne+per_O*0.00201)THEN
					Z=8
					A=18
				ELSE IF(q<per_Si+per_Mg+per_Ne+per_O*0.0024)THEN
					Z=8
					A=17
				ELSE IF(q<per_Si+per_Mg+per_Ne+per_O)THEN
					Z=8
					A=16
				ELSE IF(q<per_Si+per_Mg+per_Ne+per_O+per_N*0.00366)THEN
					Z=7
					A=15
				ELSE IF(q<per_Si+per_Mg+per_Ne+per_O+per_N)THEN
					Z=7
					A=14
				ELSE IF(q<per_Si+per_Mg+per_Ne+per_O+per_N+per_C*0.11)THEN
					Z=6
					A=13
				ELSE
					Z=6
					A=12
				END IF
			CASE('Iron')
				IF(q<(per_Ti44))THEN
					Z=22
					A=44
				ELSE IF(q<(per_Fe60+per_Ti44))THEN
					Z=26
					A=60
				ELSE IF(q<(per_Ni*(0.009256)+per_Fe60+per_Ti44))THEN
					Z=28
					A=64
				ELSE IF(q<(per_Ni*(0.009256+0.036345)+per_Fe60+per_Ti44))THEN
					Z=28
					A=62
				ELSE IF(q<(per_Ni*(0.009256+0.036345+0.011399)+per_Fe60+per_Ti44))THEN
					Z=28
					A=61
 				ELSE IF(q<(per_Ni*(0.009256+0.036345+0.011399+0.262231)+per_Fe60+per_Ti44))THEN
					Z=28
					A=60
				ELSE IF(q<(per_Ni+per_Fe60+per_Ti44))THEN
					Z=28
					A=58
 				ELSE IF(q<(per_Fe*(0.00282)+per_Ni+per_Fe60+per_Ti44))THEN
					Z=26
					A=58
 				ELSE IF(q<(per_Fe*(0.00282+0.02119)+per_Ni+per_Fe60+per_Ti44))THEN
					Z=26
					A=57
				ELSE IF(q<(per_Fe*(0.00282+0.02119+0.91754)+per_Ni+per_Fe60+per_Ti44))THEN
					Z=26
					A=56
 				ELSE
					Z=26
					A=54
				END IF
			CASE('PAMELA')
				IF(q<(per_e))THEN
					Z=-1
					A=0
 				ELSE
 					Z=1
					A=1
				END IF
		END SELECT 
	END SUBROUTINE initial_particle

	SUBROUTINE source(coords,omega)
		IMPLICIT NONE
		REAL(8),INTENT(out)::coords(3),omega
		INTEGER::arm
		REAL(8)::q,qq,qqq,b,r,phi,x,y
		CALL random_number(q)
		IF(q<per_disc_type1a)THEN
			DO
				CALL random_number(qq)
				r=-4.5*LOG(1-SQRT(qq))				!see Ferriere 2001 for SN distribution data
				IF(r<r_galaxy)EXIT
			END DO
			CALL random_number(qq)
			phi=2*pi*qq
			coords(1)=r*COS(phi)
			coords(2)=r*sin(phi)
			DO
				CALL random_number(qq)
				coords(3)=-0.325*LOG(qq)
				CALL random_number(qq)
				IF(qq<0.5)coords(3)=-coords(3)
				IF(ABS(coords(3))<z_galaxy)EXIT
			END DO
			omega=0
		ELSE IF(q<(per_disc_type1a+per_disc_type2))THEN
			DO
				DO
					CALL random_number(qq)
					r=qq*15
					CALL random_number(qq)
					y=qq*2335
					IF(r<3.7)THEN
						IF(y<602.2433*r*EXP(-((r-3.7)/2.1)**2))EXIT
					ELSE
						IF(y<809.3358*r*EXP(-(r/6.8)**2))EXIT
					END IF
				END DO
				IF(r<r_galaxy)EXIT
			END DO
			CALL random_number(qq)
			phi=2*pi*qq
			coords(1)=r*COS(phi)
			coords(2)=r*sin(phi)
			DO
				CALL random_number(qq)
				IF(qq<0.79)THEN
					CALL gasdev(qqq)
					coords(3)=0.212*qqq
				ELSE
					CALL gasdev(qqq)
					coords(3)=0.636*qqq
				END IF
				IF(ABS(coords(3))<z_galaxy)EXIT
			END DO
 			omega=0
		ELSE IF(q<(per_disc_type1a+per_disc_type2+per_2arm_set))THEN
			DO
				CALL gasdev(qq)
				r=6.8*ABS(qq)
				IF(r>inner_r_2arm_set.and.r<r_galaxy)EXIT
			END DO
			b=1/TAN(pi/2-pitch_angle_2arm_set)
			phi=LOG(r/inner_r_2arm_set)/b
			CALL random_number(qq)
			arm=(qq*2)/1
			phi=phi+pi*arm+start_angle_2arm_set
			coords(1)=r*COS(phi)
			coords(2)=r*SIN(phi)
			DO
				CALL random_number(qq)
				IF(qq<0.79)THEN
					CALL gasdev(qqq)
					coords(3)=0.212*qqq
				ELSE
					CALL gasdev(qqq)
					coords(3)=0.636*qqq
				END IF
				IF(ABS(coords(3))<z_galaxy)EXIT
			END DO
			CALL random_number(qq)
			x=width_2arm_set*(2*qq-1)
			coords(1)=x*SIN(phi)+coords(1)
			coords(2)=x*COS(phi)+coords(2)
			omega=omega_2arm_set
		ELSE IF(q<(per_disc_type1a+per_disc_type2+per_2arm_set+per_4arm_set))THEN
			DO
				CALL gasdev(qq)
				r=6.8*ABS(qq)
				IF(r>inner_r_4arm_set.and.r<r_galaxy)EXIT
			END DO
			b=1/TAN(pi/2-pitch_angle_4arm_set)
			phi=LOG(r/inner_r_4arm_set)/b
			CALL random_number(qq)
			arm=(qq*4)/1
			phi=phi+pi*arm/2+start_angle_4arm_set
			coords(1)=r*COS(phi)
			coords(2)=r*SIN(phi)
			DO																		
				CALL random_number(qq)
				IF(qq<0.79)THEN
					CALL gasdev(qqq)
					coords(3)=0.212*qqq
				ELSE
					CALL gasdev(qqq)
					coords(3)=0.636*qqq
				END IF
				IF(ABS(coords(3))<z_galaxy)EXIT
			END DO
			CALL random_number(qq)
			x=width_4arm_set*(2*qq-1)
			coords(1)=x*SIN(phi)+coords(1)
			coords(2)=x*COS(phi)+coords(2)
			omega=omega_4arm_set
		END IF
	END SUBROUTINE source

	SUBROUTINE step(dt,l,omega,disc_r,disc_v,coords)
		IMPLICIT NONE
		REAL(8),INTENT(in)::dt,l,omega,disc_r(1:n),disc_v(1:n)
		REAL(8),INTENT(in out)::coords(3)
		REAL(8)::r,v,dt_sec,q,g,teta,phi
		r=SQRT(coords(1)**2+coords(2)**2)
		phi=ATAN(coords(2)/coords(1))
		IF(coords(1)<0)THEN
			phi=phi+pi
		END IF
		v=interp(n,disc_r,disc_v,r)
		dt_sec=dt*megayear_sec
		phi=phi-(v/(r*kpc_km)-omega/kpc_km)*dt_sec
		CALL gasdev(q)
		coords(1)=r*COS(phi)+l*q
		CALL gasdev(q)
		coords(2)=r*SIN(phi)+l*q
		CALL gasdev(q)
		coords(3)=coords(3)+l*q
	END SUBROUTINE step
END MODULE diffusion