MODULE spalation
	USE input
	USE interpol
	USE gauss
	IMPLICIT NONE
	INTEGER::A,Z

	REAL(8)::half_t_Ni59_rel,half_t_Fe60_rel,half_t_Mn53_rel,half_t_Al26_rel,half_t_Be10_rel,alpha,beta
	CONTAINS

	SUBROUTINE particle_density(e,initial_energy,energy,v,g,coords,weight,nuc_density,events,dt,l)
		IMPLICIT NONE

		INTEGER,INTENT(in)::e
		REAL(8),INTENT(in)::initial_energy,energy,v,g,coords(3),weight

		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(8),INTENT(out)::dt,l

		REAL(8)::p,D,t

		IF(A==0)THEN
			D=D0*v*(energy/3)**(delta)
		ELSE
			p=A*SQRT(energy**2+2*energy)
			D=D0*v*(p/(3*Z))**(delta)
		END IF
		t=z_galaxy**2/(2*D)
		dt=0.01*t
		l=SQRT(2*D*dt)

		SELECT CASE(model)
			CASE('disc')
				IF(SQRT((SQRT(coords(1)**2+coords(2)**2)-8.5)**2+coords(3)**2)<0.15)THEN
					dt=0.1*dt
					l=SQRT(2*D*dt)
					IF(SQRT((SQRT(coords(1)**2+coords(2)**2)-8.5)**2+coords(3)**2)<0.12)THEN
						dt=0.1*dt
						l=SQRT(2*D*dt)
						IF(SQRT((SQRT(coords(1)**2+coords(2)**2)-8.5)**2+coords(3)**2)<0.1)THEN
							dt=0.16*dt
							l=SQRT(2*D*dt)
						END IF
					END IF
				END IF
				IF(SQRT((SQRT(coords(1)**2+coords(2)**2)-8.5)**2+coords(3)**2)<0.1)THEN
					nuc_density(Z,A,e)=nuc_density(Z,A,e)+dt*(alpha*initial_energy**(-(beta-delta)))*weight		!observations see dN/dE~-2.7, my beans are logarithmic->dN/d(log(E))=E*dN/dE~-1.7, need to subtract the delta which comes from diffusion
					events(Z,A,e)=events(Z,A,e)+1
				END IF
			CASE('spiral arms')
				IF(SQRT((ABS(coords(1))-8.5)**2+coords(2)**2+coords(3)**2)<0.15)THEN
					dt=0.1*dt
					l=SQRT(2*D*dt)
					IF(SQRT((ABS(coords(1))-8.5)**2+coords(2)**2+coords(3)**2)<0.12)THEN
						dt=0.1*dt
						l=SQRT(2*D*dt)
						IF(SQRT((ABS(coords(1))-8.5)**2+coords(2)**2+coords(3)**2)<0.1)THEN
							dt=0.16*dt
							l=SQRT(2*D*dt)
						END IF
					END IF
				END IF
				IF(SQRT((ABS(coords(1))-8.5)**2+coords(2)**2+coords(3)**2)<0.1)THEN
					nuc_density(Z,A,e)=nuc_density(Z,A,e)+dt*(alpha*initial_energy**(-(beta-delta)))*weight
					events(Z,A,e)=events(Z,A,e)+1
				END IF
		END SELECT 

	END SUBROUTINE particle_density
	
	SUBROUTINE break_particle(e,g,energy,xsec,e_sec,break)
		IMPLICIT NONE

		INTEGER,INTENT(in)::e
	
		REAL(8),INTENT(in)::g,energy
		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)

		INTEGER,INTENT(out)::e_sec
		LOGICAL,INTENT(out)::break

		INTEGER::jz,ja
		REAL(8)::q,tot_xsec
		LOGICAL::particle_broke

		CALL random_number(q)

		SELECT CASE(plot)
			CASE('B/C','Iron')
				IF(q<g*xsec(Z,A,e,min_Z,max_A,e))break=.true.

				particle_broke=.false.
				tot_xsec=xsec(Z,A,e,min_Z,max_A,e)

				DO jz=min_Z,Z

					DO ja=min_A,A-1

						tot_xsec=tot_xsec+xsec(Z,A,e,jz,ja,e)

						IF(q<g*tot_xsec)THEN
							Z=jz
							A=ja

							particle_broke=.true.
							EXIT

						END IF
	
					END DO

					IF(particle_broke)EXIT
				END DO

			CASE('PAMELA')
				particle_broke=.false.
				tot_xsec=0								!*******need to add breaking below min_E*******
				IF(Z==1.AND.A==1)THEN
					DO ja=min_A,max_A
						DO jz=min_Z,max_Z
							DO e_sec=min_cell_E,e		!*******need to add min energy*****
								tot_xsec=tot_xsec+xsec(Z,A,e,jz,ja,e_sec)

								IF(q<g*tot_xsec)THEN
									Z=jz
									A=ja

									break=.true.

									particle_broke=.true.
									EXIT

								END IF

							END DO

							IF(particle_broke)EXIT

						END DO
						IF(particle_broke)EXIT

					END DO
				END IF
		END SELECT
	END SUBROUTINE break_particle

	
	SUBROUTINE particle_spalation(e,initial_energy,energy,v,g,coords,weight,xsec,nuc_density,events,e_sec,dt,l,break)
		IMPLICIT NONE
		INTEGER,INTENT(in)::e
		REAL(8),INTENT(in)::initial_energy,energy,v,g,coords(3),weight
		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)
		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)
		INTEGER,INTENT(out)::e_sec
		REAL(8),INTENT(out)::dt,l
		LOGICAL,INTENT(out)::break
		REAL(8)::q,pEC

		CALL random_number(q)
		pEC=(ECa*(Z**(4.5)))/(((energy/0.5))**(ECdelta))
		IF(Z==4.and.A==10)THEN

			IF(q<(1-EXP(-dt/half_t_Be10_rel)))Z=5
		ELSE IF(Z==13.and.A==26)THEN

			IF(q<(1-EXP(-dt/half_t_Al26_rel)))Z=12			


		ELSE IF(Z==26.and.A==60)THEN

			IF(q<(1-EXP(-dt/half_t_Fe60_rel)))Z=27
		ELSE IF(Z==22.and.A==44)THEN

			IF(q<pEC*g)Z=20

		ELSE IF(Z==23.and.A==49)THEN

			IF(q<pEC*g)Z=22

		ELSE IF(Z==24.and.A==51)THEN

			IF(q<pEC*g)Z=23

		ELSE IF(Z==25.and.A==54)THEN

			IF(q<pEC*g)THEN
				Z=24
			ELSE
				Z=26
			END IF

		ELSE IF(Z==26.and.A==55)THEN

			IF(q<pEC*g)Z=25

		ELSE IF(Z==27.and.A==57)THEN

			IF(q<pEC*g)Z=26

		END IF

		CALL particle_density(e,initial_energy,energy,v,g,coords,weight,nuc_density,events,dt,l)

		CALL break_particle(e,g,energy,xsec,e_sec,break)
	END SUBROUTINE particle_spalation
END MODULE spalation
