MODULE output
	USE input
	USE interpol
	IMPLICIT NONE
	CONTAINS
		SUBROUTINE plots(nuc_density,density)
		IMPLICIT NONE
		REAL(8),INTENT(in)::nuc_density(min_Z:max_Z,min_A:max_A,min_cell_E:max_cell_E)
		REAL(4),INTENT(in)::density(-cell_XY:cell_XY,-cell_XY:cell_XY,-cell_Z:cell_Z,min_cell_E:max_cell_E)
		SELECT CASE(calculate)
		
			CASE('S/P')
		
				CALL plot_nuc(nuc_density)
			CASE('galaxy')
				SELECT CASE(density_plot)
					CASE('gamma')
						CALL plot_gamma(density)
					CASE('X-Y')
						CALL plot_xy(density)
					CASE('X-Z')
						CALL plot_xz(density)
					CASE('Z-phi')
						CALL plot_z_phi(density)
			
					CASE('anisotropy')

						CALL plot_anisotropy(density)
				END SELECT 
		END SELECT 
	END SUBROUTINE plots
	
	SUBROUTINE plot_nuc(nuc_density)
		IMPLICIT NONE
		REAL(8),INTENT(in)::nuc_density(min_Z:max_Z,min_A:max_A,min_cell_E:max_cell_E)

		REAL(8)::O,N,C,B,Be,Be10,Be9

		REAL(8)::Sc,Ti,V,Cr,Mn,Fe,Co,Ni,Ca44,Ti44,Ti49,V49,V51,Cr51,Fe57,Co57,Fe56,Fe60
		REAL(8)::electrons,positrons,protons
		INTEGER::e,ia

		SELECT CASE(plot)

			CASE('B/C')
				OPEN(UNIT=100,FILE='BC.out')
				OPEN(UNIT=101,FILE='CO.out')
				OPEN(UNIT=102,FILE='NO.out')

				OPEN(UNIT=103,FILE='Be10Be9.out')
				OPEN(UNIT=104,FILE='Be.out')
				OPEN(UNIT=105,FILE='B.out')
				OPEN(UNIT=106,FILE='C.out')
				OPEN(UNIT=107,FILE='N.out')
				OPEN(UNIT=108,FILE='O.out')
				DO e=min_cell_E,max_cell_E
					O=0
	
					N=0
	
					C=0

					B=0
					Be=0

					Be10=0

					Be9=0

					DO ia=min_A,max_A

						O=O+nuc_density(8,ia,e)
						N=N+nuc_density(7,ia,e)
						C=C+nuc_density(6,ia,e)

						B=B+nuc_density(5,ia,e)

						Be=Be+nuc_density(4,ia,e)

					END DO

					Be10=nuc_density(4,10,e)

					Be9=nuc_density(4,9,e)
					WRITE(100,*)B/C
					WRITE(101,*)C/O
					WRITE(102,*)N/O
					WRITE(103,*)Be10/Be9
					WRITE(104,*)Be
					WRITE(105,*)B
					WRITE(106,*)C
					WRITE(107,*)N
					WRITE(108,*)O
				END DO
				CLOSE(100)
				CLOSE(101)
				CLOSE(102)
				CLOSE(103)
				CLOSE(104)
				CLOSE(105)
				CLOSE(106)
				CLOSE(107)
				CLOSE(108)
			CASE('Iron')

				OPEN(UNIT=100,FILE='SubIron_Iron3.out')
				OPEN(UNIT=101,FILE='Ni_Fe.out')

				OPEN(UNIT=102,FILE='Ca44Ti44.out')

				OPEN(UNIT=103,FILE='Ti44Ti.out')

				OPEN(UNIT=104,FILE='Ti49V49.out')
				OPEN(UNIT=105,FILE='V51Cr51.out')

				OPEN(UNIT=106,FILE='Fe57Co57.out')

				OPEN(UNIT=107,FILE='Fe60Fe56.out')

				OPEN(UNIT=108,FILE='Sc.out')

				OPEN(UNIT=109,FILE='Ti.out')

				OPEN(UNIT=110,FILE='V.out')

				OPEN(UNIT=111,FILE='Cr.out')
				OPEN(UNIT=112,FILE='Mn.out')
				OPEN(UNIT=113,FILE='Fe.out')
				OPEN(UNIT=114,FILE='Co.out')
				OPEN(UNIT=115,FILE='Ni.out')

				OPEN(UNIT=116,FILE='Ca44.out')

				OPEN(UNIT=117,FILE='Ti44.out')

				OPEN(UNIT=118,FILE='Ti49.out')

				OPEN(UNIT=119,FILE='V49.out')

				OPEN(UNIT=120,FILE='V51.out')

				OPEN(UNIT=121,FILE='Cr51.out')

				OPEN(UNIT=122,FILE='Fe56.out')

				OPEN(UNIT=123,FILE='Fe60.out')

				DO e=min_cell_E,max_cell_E
					Ni=0

					Co=0

					Fe=0

					Mn=0

					Cr=0

					V=0

					Ti=0

					Sc=0

					DO ia=min_A,max_A
						Ni=Ni+nuc_density(28,ia,e)

						Co=Co+nuc_density(27,ia,e)

						Fe=Fe+nuc_density(26,ia,e)

						Mn=Mn+nuc_density(25,ia,e)

						Cr=Cr+nuc_density(24,ia,e)

						V=V+nuc_density(23,ia,e)

						Ti=Ti+nuc_density(22,ia,e)

						Sc=Sc+nuc_density(21,ia,e)

					END DO

					Ca44=nuc_density(20,44,e)

					Ti44=nuc_density(22,44,e)

					Ti49=nuc_density(22,49,e)
					V49=nuc_density(23,49,e)

					V51=nuc_density(23,51,e)

					Cr51=nuc_density(24,51,e)

					Fe57=nuc_density(26,57,e)

					Co57=nuc_density(27,57,e)

					Fe56=nuc_density(26,56,e)

					Fe60=nuc_density(26,60,e)

					WRITE(100,*)(V+Ti+Sc)/Fe
					WRITE(101,*)Ni/Fe
					WRITE(102,*)Ca44/Ti44
					WRITE(103,*)Ti44/Ti
					WRITE(104,*)Ti49/V49

					WRITE(105,*)V51/Cr51
					WRITE(106,*)Fe57/Co57
					WRITE(107,*)Fe60/Fe56
					WRITE(108,*)Sc
					WRITE(109,*)Ti
					WRITE(110,*)V
					WRITE(111,*)Cr
					WRITE(112,*)Mn
					WRITE(113,*)Fe
					WRITE(114,*)Co
					WRITE(115,*)Ni
					WRITE(116,*)Ca44
					WRITE(117,*)Ti44
					WRITE(118,*)Ti49
					WRITE(119,*)V49
					WRITE(120,*)V51
					WRITE(121,*)Cr51
					WRITE(122,*)Fe56
					WRITE(123,*)Fe60
				END DO

				CLOSE(100)
				CLOSE(101)
				CLOSE(102)
				CLOSE(103)
				CLOSE(104)
				CLOSE(105)
				CLOSE(106)
				CLOSE(107)
				CLOSE(108)
				CLOSE(109)
				CLOSE(110)
				CLOSE(111)
				CLOSE(112)
				CLOSE(113)
				CLOSE(114)
				CLOSE(115)
				CLOSE(116)
				CLOSE(117)
				CLOSE(118)
				CLOSE(119)
				CLOSE(120)
				CLOSE(121)
				CLOSE(122)
				CLOSE(123)
			CASE('PAMELA')
				OPEN(UNIT=100,FILE='positron fraction.out')
				OPEN(UNIT=101,FILE='electrons.out')
				OPEN(UNIT=102,FILE='positrons.out')

				OPEN(UNIT=103,FILE='protons.out')
				OPEN(UNIT=104,FILE='EP.out')
				DO e=min_cell_E,max_cell_E
					electrons=0

					positrons=0

					protons=0

					positrons=nuc_density(1,0,e)
					electrons=nuc_density(-1,0,e)
					protons=nuc_density(1,1,e)

					WRITE(100,*)positrons/(positrons+electrons)
					WRITE(101,*)electrons
					WRITE(102,*)positrons
					WRITE(103,*)protons
					WRITE(104,*)electrons/protons
				END DO
				CLOSE(100)
				CLOSE(101)
				CLOSE(102)
				CLOSE(103)
				CLOSE(104)
		END SELECT 
	END SUBROUTINE plot_nuc
	
	SUBROUTINE calculate_chai(nuc_density,events)

		IMPLICIT NONE
		INTEGER,INTENT(in)::events(min_Z:max_Z,min_A:max_A,min_cell_E:max_cell_E)
		REAL(8),INTENT(in)::nuc_density(min_Z:max_Z,min_A:max_A,min_cell_E:max_cell_E)
		INTEGER::A,Z,i,j,m,num_points
		REAL(8)::energy_obs,data_obs,error_obs,data_data,data_data1,data_data2,error_data,error_data1,error_data2
		REAL(8)::chai,null,obser(1:12)
		CHARACTER::data_name*15
		REAL(8),ALLOCATABLE::data_sim(:),energy_array(:),nuc_Z(:,:),num(:,:),num_sqr(:,:),Var(:,:)
		REAL(8),ALLOCATABLE::EC1_data_sim(:),EC2_data_sim(:),energy_array_solar(:)
		REAL(8),ALLOCATABLE::error(:,:,:),error_Z(:,:),error_sub(:),EC1_error_44(:),EC2_error_44(:)
		REAL(8),ALLOCATABLE::EC1_error_49(:),EC2_error_49(:),EC1_error_51(:),EC2_error_51(:)
		REAL(8),ALLOCATABLE::EC1_error_60(:),EC2_error_60(:)
		SELECT CASE(plot)

			CASE('B/C')
				OPEN(UNIT=100,FILE='BC.out')
				OPEN(UNIT=700,FILE='BoverC.csv')
			CASE('Iron')

				OPEN(UNIT=100,FILE='SubIron_Iron3.out')
				OPEN(UNIT=101,FILE='Ca44.out')
				OPEN(UNIT=102,FILE='Ti44.out')
				OPEN(UNIT=103,FILE='Ti49.out')
				OPEN(UNIT=104,FILE='V49.out')
				OPEN(UNIT=105,FILE='V51.out')
				OPEN(UNIT=106,FILE='Cr51.out')
				OPEN(UNIT=107,FILE='Fe60.out')
				OPEN(UNIT=108,FILE='Fe56.out')
				OPEN(UNIT=700,FILE='subIron.csv')
				OPEN(UNIT=701,FILE='44.csv')
				OPEN(UNIT=702,FILE='49.csv')
				OPEN(UNIT=703,FILE='51.csv')
				OPEN(UNIT=704,FILE='60.csv')
		END SELECT
		OPEN(UNIT=10,FILE='energy.out')
		i=max_cell_E-min_cell_E+1
		ALLOCATE(energy_array(1:i))
		energy_array=0
		ALLOCATE(energy_array_solar(1:i))
		energy_array_solar=0
		ALLOCATE(data_sim(1:i))
		data_sim=0
		ALLOCATE(EC1_data_sim(1:i))
		EC1_data_sim=0
		ALLOCATE(EC2_data_sim(1:i))
		EC2_data_sim=0
		ALLOCATE(nuc_Z(min_Z:max_Z,1:i))
		nuc_Z=0
		ALLOCATE(num(min_Z:max_Z,min_A:max_A))
		num=0
		ALLOCATE(num_sqr(min_Z:max_Z,min_A:max_A))
		num_sqr=0
		ALLOCATE(Var(min_Z:max_Z,min_A:max_A))
		Var=0
		ALLOCATE(error(min_Z:max_Z,min_A:max_A,1:i))
		error=0
		ALLOCATE(error_Z(min_Z:max_Z,1:i))
		error_Z=0
		ALLOCATE(error_sub(1:i))
		error_sub=0
		ALLOCATE(EC1_error_44(1:i))
		EC1_error_44=0
		ALLOCATE(EC2_error_44(1:i))
		EC2_error_44=0
		ALLOCATE(EC1_error_49(1:i))
		EC1_error_49=0
		ALLOCATE(EC2_error_49(1:i))
		EC2_error_49=0
		ALLOCATE(EC1_error_51(1:i))
		EC1_error_51=0
		ALLOCATE(EC2_error_51(1:i))
		EC2_error_51=0
		ALLOCATE(EC1_error_60(1:i))
		EC1_error_60=0
		ALLOCATE(EC2_error_60(1:i))
		EC2_error_60=0
		A=0
		Z=0
		DO j=1,i
			READ(10,*)energy_array(j)
			READ(100,*)data_sim(j)
			Do Z=min_Z,max_Z
				Do A=min_A,max_A
					nuc_Z(Z,j)=nuc_Z(Z,j)+nuc_density(Z,A,j+min_cell_E-1)
					num(Z,A)=num(Z,A)+nuc_density(Z,A,j+min_cell_E-1)
					num_sqr(Z,A)=num_sqr(Z,A)+nuc_density(Z,A,j+min_cell_E-1)**(2.)
				END DO
			END DO
		END DO
		CLOSE(10)
		CLOSE(100)
		num=num/i
		num_sqr=num_sqr/i
		Do Z=min_Z,max_Z
			Do A=min_A,max_A
				Var(Z,A)=SQRT(num_sqr(Z,A)-num(Z,A)**(2.))
			END DO
		END DO
		DO j=1,i
			Do Z=min_Z,max_Z
				Do A=min_A,max_A
					IF(events(Z,A,j+min_cell_E-1)>0)THEN
						error(Z,A,j)=Var(Z,A)/SQRT(REAL(events(Z,A,j+min_cell_E-1)))
						error_Z(Z,j)=error_Z(Z,j)+error(Z,A,j)
					END IF
				END DO
			END DO
		END DO
		SELECT CASE(plot)

			CASE('B/C')
				DO j=1,i
					error_sub(j)=(error_Z(5,j)/nuc_Z(6,j))+((nuc_Z(5,j)*error_Z(6,j))/(nuc_Z(6,j)**(2.)))
				END DO
			CASE('Iron')

				DO j=1,i
					error_sub(j)=((error_Z(21,j)+error_Z(22,j)+error_Z(23,j))/nuc_Z(26,j))+&
							&(((nuc_Z(21,j)+nuc_Z(22,j)+nuc_Z(23,j))*error_Z(26,j))/(nuc_Z(26,j)**(2.)))
					EC1_error_44(j)=error(20,44,j)
					EC2_error_44(j)=error(22,44,j)
					EC1_error_49(j)=error(22,49,j)
					EC2_error_49(j)=error(23,49,j)
					EC1_error_51(j)=error(23,51,j)
					EC2_error_51(j)=error(24,51,j)
					EC1_error_51(j)=error(26,60,j)
					EC2_error_51(j)=error(26,56,j)
				END DO
		END SELECT
		m=0
		READ(700,*)			!first line is head
		DO
			READ(700,*,END=900)
			m=m+1
		END DO
		900 REWIND(700)
		chai=0
		num_points=0
		READ(700,*)
		DO j=1,m
			READ(700,*)data_name,obser(1),obser(2),obser(3),obser(4),obser(5),obser(6),obser(7),obser(8),obser(9),&
						&obser(10),obser(11),obser(12)
			SELECT CASE(plot)
				CASE('B/C')	
					SELECT CASE(data_name)
						CASE('AMS-02')
							IF(obser(1)>=energy_start.AND.obser(1)<=energy_end)THEN
								num_points=num_points+1
								error_data=interp(i,energy_array,error_sub,obser(1))
								data_data=interp(i,energy_array,data_sim,obser(1))
								chai=chai+((obser(5)-data_data)**2)/((obser(6)+error_data)**2)
							END IF
					END SELECT
				CASE('Iron')	
					SELECT CASE(data_name)
						CASE('SANRIKU','HEAO-3')
							IF(obser(1)>=energy_start.AND.obser(1)<=energy_end)THEN
								num_points=num_points+1
								error_data=interp(i,energy_array,error_sub,obser(1))
								data_data=interp(i,energy_array,data_sim,obser(1))
								chai=chai+((obser(5)-data_data)**2)/((obser(6)+error_data)**2)
							END IF
					END SELECT
			END SELECT
		END DO
		chai=chai/num_points
		write(*,*)plot,'chi square=',chai
		CLOSE(700)
		SELECT CASE(plot)
			CASE('Iron')
				DO j=1,i
					READ(101,*)EC1_data_sim(j)
					READ(102,*)EC2_data_sim(j)
				END DO
				CLOSE(101)
				CLOSE(102)
				m=0
				READ(701,*)			!first line is head
				DO
					READ(701,*,END=901)
					m=m+1
				END DO
				901 REWIND(701)
				chai=0
				num_points=0
				READ(701,*)
				DO j=1,m
					READ(701,*)data_name,obser(1),obser(2),obser(3),obser(4),obser(5),obser(6),obser(7),obser(8),obser(9),&
						&obser(10),obser(11),obser(12)
					SELECT CASE(data_name)
						CASE('CRIS-solar_min','CRIS-solar_max')
							IF(obser(1)>=energy_start.AND.obser(1)<=energy_end)THEN
								num_points=num_points+1
								energy_array_solar=energy_array-(20./44.)*(obser(12)/1000)		!solar modulation in MV
								error_data1=interp(i,energy_array_solar,EC1_error_44,obser(1))
								data_data1=interp(i,energy_array_solar,EC1_data_sim,obser(1))
								energy_array_solar=energy_array-(22./44.)*(obser(12)/1000)		!solar modulation in MV
								error_data2=interp(i,energy_array_solar,EC2_error_44,obser(1))
								data_data2=interp(i,energy_array_solar,EC2_data_sim,obser(1))
								error_data=(error_data1/data_data2)+((error_data2*data_data1)/(data_data2**(2.)))
								data_data=data_data1/data_data2
								chai=chai+((obser(5)-data_data)**2)/((obser(6)+error_data)**2)
							END IF
					END SELECT
				END DO
				CLOSE(701)
				write(*,*)'EC44             ','chi square=',chai,'number of data points=',num_points
				
				DO j=1,i
					READ(103,*)EC1_data_sim(j)
					READ(104,*)EC2_data_sim(j)
				END DO
				CLOSE(103)
				CLOSE(104)
				m=0
				READ(702,*)			!first line is head
				DO
					READ(702,*,END=902)
					m=m+1
				END DO
				902 REWIND(702)
				chai=0
				num_points=0
				READ(702,*)
				DO j=1,m
					READ(702,*)data_name,obser(1),obser(2),obser(3),obser(4),obser(5),obser(6),obser(7),obser(8),obser(9),&
						&obser(10),obser(11),obser(12)
					SELECT CASE(data_name)
						CASE('CRIS-solar_min','CRIS-solar_max')
							IF(obser(1)>=energy_start.AND.obser(1)<=energy_end)THEN
								num_points=num_points+1
								energy_array_solar=energy_array-(22./49.)*(obser(12)/1000)		!solar modulation in MV
								error_data1=interp(i,energy_array_solar,EC1_error_49,obser(1))
								data_data1=interp(i,energy_array_solar,EC1_data_sim,obser(1))
								energy_array_solar=energy_array-(23./49.)*(obser(12)/1000)		!solar modulation in MV
								error_data2=interp(i,energy_array_solar,EC2_error_49,obser(1))
								data_data2=interp(i,energy_array_solar,EC2_data_sim,obser(1))
								error_data=(error_data1/data_data2)+((error_data2*data_data1)/(data_data2**(2.)))
								data_data=data_data1/data_data2
								chai=chai+((obser(5)-data_data)**2)/((obser(6)+error_data)**2)
							END IF
					END SELECT
				END DO
				CLOSE(702)
				write(*,*)'EC49             ','chi square=',chai,'number of data points=',num_points
				
				DO j=1,i
					READ(105,*)EC1_data_sim(j)
					READ(106,*)EC2_data_sim(j)
				END DO
				CLOSE(105)
				CLOSE(106)
				m=0
				READ(703,*)			!first line is head
				DO
					READ(703,*,END=903)
					m=m+1
				END DO
				903 REWIND(703)
				chai=0
				num_points=0
				READ(703,*)
				DO j=1,m
					READ(703,*)data_name,obser(1),obser(2),obser(3),obser(4),obser(5),obser(6),obser(7),obser(8),obser(9),&
						&obser(10),obser(11),obser(12)
					SELECT CASE(data_name)
						CASE('CRIS-solar_min','CRIS-solar_max')
							IF(obser(1)>=energy_start.AND.obser(1)<=energy_end)THEN
								num_points=num_points+1
								energy_array_solar=energy_array-(23./51.)*(obser(12)/1000)		!solar modulation in MV
								error_data1=interp(i,energy_array_solar,EC1_error_51,obser(1))
								data_data1=interp(i,energy_array_solar,EC1_data_sim,obser(1))
								energy_array_solar=energy_array-(24./51.)*(obser(12)/1000)		!solar modulation in MV
								error_data2=interp(i,energy_array_solar,EC2_error_51,obser(1))
								data_data2=interp(i,energy_array_solar,EC2_data_sim,obser(1))
								error_data=(error_data1/data_data2)+((error_data2*data_data1)/(data_data2**(2.)))
								data_data=data_data1/data_data2
								chai=chai+((obser(5)-data_data)**2)/((obser(6)+error_data)**2)
							END IF
					END SELECT
				END DO
				chai=chai
				write(*,*)'EC51             ','chi square=',chai,'number of data points=',num_points
				CLOSE(703)
				
				DO j=1,i
					READ(107,*)EC1_data_sim(j)
					READ(108,*)EC2_data_sim(j)
				END DO
				CLOSE(107)
				CLOSE(108)
				m=0
				READ(704,*)			!first line is head
				DO
					READ(704,*,END=904)
					m=m+1
				END DO
				904 REWIND(704)
				chai=0
				num_points=0
				READ(704,*)
				DO j=1,m
					READ(704,*)data_name,obser(1),obser(2),obser(3),obser(4),obser(5),obser(6),obser(7),obser(8),obser(9),&
						&obser(10),obser(11),obser(12)
					SELECT CASE(data_name)
						CASE('CRIS')
							IF(obser(1)>=energy_start.AND.obser(1)<=energy_end)THEN
								num_points=num_points+1
								energy_array_solar=energy_array-(26./60.)*(obser(12)/1000)		!solar modulation in MV
								error_data1=interp(i,energy_array_solar,EC1_error_60,obser(1))
								data_data1=interp(i,energy_array_solar,EC1_data_sim,obser(1))
								energy_array_solar=energy_array-(26./56.)*(obser(12)/1000)		!solar modulation in MV
								error_data2=interp(i,energy_array_solar,EC2_error_60,obser(1))
								data_data2=interp(i,energy_array_solar,EC2_data_sim,obser(1))
								error_data=(error_data1/data_data2)+((error_data2*data_data1)/(data_data2**(2.)))
								data_data=data_data1/data_data2
								chai=chai+((obser(5)-data_data)**2)/((obser(6)+error_data)**2)
							END IF
					END SELECT
				END DO
				chai=chai
				write(*,*)'EC60             ','chi square=',chai,'number of data points=',num_points
				CLOSE(704)
		END SELECT
	END SUBROUTINE calculate_chai

	SUBROUTINE plot_xy(density)
		IMPLICIT NONE
		REAL(4),INTENT(in)::density(-cell_XY:cell_XY,-cell_XY:cell_XY,-cell_Z:cell_Z,min_cell_E:max_cell_E)
		INTEGER::i,j
		REAL(8)::x,y,z
		OPEN(UNIT=20,file='density_xy.out')
		OPEN(UNIT=21,file='x_xy.out')
		OPEN(UNIT=22,file='y_xy.out')
		DO i=-cell_XY,cell_XY
			y=i*grid
			DO j=-cell_XY,cell_XY
				x=j*grid
				WRITE(20,'(1es12.4)',ADVANCE='NO')density(j,i,0,min_cell_E)
				WRITE(21,'(f8.2)',ADVANCE='NO')x
				WRITE(22,'(f8.2)',ADVANCE='NO')y
			END DO
			WRITE(20,*)""
			WRITE(21,*)""
			WRITE(22,*)""
		END DO
		CLOSE(20)
		CLOSE(21)
		CLOSE(22)
	END SUBROUTINE plot_xy

	SUBROUTINE plot_gamma(density)
		IMPLICIT NONE
		REAL(4),INTENT(in)::density(-cell_XY:cell_XY,-cell_XY:cell_XY,-cell_Z:cell_Z,min_cell_E:max_cell_E)
		INTEGER::i,j,z
		OPEN(UNIT=17,file='density.out')
		OPEN(UNIT=18,file='x_xy.out')
		OPEN(UNIT=19,file='y_xy.out')
		DO i=-cell_XY,cell_XY
			DO j=-cell_XY,cell_XY
				WRITE(18,'(f8.2)',ADVANCE='NO')j*grid
				WRITE(19,'(f8.2)',ADVANCE='NO')i*grid
			END DO
			WRITE(18,*)""
			WRITE(19,*)""
		END DO
		DO z=-cell_Z,cell_Z
			WRITE(17,*)""
			DO i=-cell_XY,cell_XY
				DO j=-cell_XY,cell_XY
					WRITE(17,'(1es12.4)',ADVANCE='NO')density(j,i,z,min_cell_E)
				END DO
				WRITE(17,*)""
			END DO
		END DO
		CLOSE(17)
		CLOSE(18)
		CLOSE(19)
	END SUBROUTINE plot_gamma

	SUBROUTINE plot_z_phi(density)
		IMPLICIT NONE
		REAL(4),INTENT(in)::density(-cell_XY:cell_XY,-cell_XY:cell_XY,-cell_Z:cell_Z,min_cell_E:max_cell_E)
		INTEGER::i,j,n,m
		REAL(8)::x,y,z,phi
		OPEN(UNIT=30,file='density_z_phi.out')
		OPEN(UNIT=31,file='z_z_phi.out')
		OPEN(UNIT=32,file='phi_z_phi.out')
		DO i=-cell_Z,cell_Z
			z=i*grid
			DO j=0,628
				phi=REAL(j)/100
				x=SQRT(r_plot**2/(1+(TAN(phi))**2))
				n=x/grid
				y=x*TAN(phi)
				m=y/grid
				WRITE(30,'(1es12.4)',ADVANCE='NO')density(n,m,i,min_cell_E)
				WRITE(31,'(f8.2)',ADVANCE='NO')z
				WRITE(32,'(f8.2)',ADVANCE='NO')phi
			END DO
			WRITE(30,*)""
			WRITE(31,*)""
			WRITE(32,*)""
		END DO
		CLOSE(30)
		CLOSE(31)
		CLOSE(32)
	END SUBROUTINE plot_z_phi

	SUBROUTINE plot_xz(density)
		IMPLICIT NONE
		REAL(4),INTENT(in)::density(-cell_XY:cell_XY,-cell_XY:cell_XY,-cell_Z:cell_Z,min_cell_E:max_cell_E)
		INTEGER::i,j
		REAL(8)::x,y,z
		OPEN(UNIT=40,file='density_xz.out')
		OPEN(UNIT=41,file='x_xz.out')
		OPEN(UNIT=42,file='y_xz.out')
		DO i=-cell_Z,cell_Z
			z=i*grid
			DO j=-cell_XY,cell_XY
				x=j*grid
				WRITE(40,'(1es12.4)',ADVANCE='NO')density(j,0,i,min_cell_E)
				WRITE(41,'(f8.2)',ADVANCE='NO')x
				WRITE(42,'(f8.2)',ADVANCE='NO')z
			END DO
			WRITE(40,*)""
			WRITE(41,*)""
			WRITE(42,*)""
		END DO
		CLOSE(40)
		CLOSE(41)
		CLOSE(42)
	END SUBROUTINE plot_xz


	
	SUBROUTINE plot_anisotropy(density)

		IMPLICIT NONE
		REAL(4), INTENT(in)::density(-cell_XY:cell_XY,-cell_XY:cell_XY,-cell_Z:cell_Z,min_cell_E:max_cell_E)
		INTEGER::e,x,y,z,x_min,x_max,y_min,y_max,z_min,z_max

		REAL(8)::grad,gradx,grady,gradz
		x=8.5/grid
		y=0/grid
		z=0/grid

		x_min=x-1

		x_max=x+1

		y_min=y-1

		y_max=y+1

		z_min=z-1

		z_max=z+1

		OPEN(UNIT=50,file='anisotropy.out')

		WRITE(50,'(6A40)')'energy[GeV]','gradient X direction[n/kpc]','gradient Y direction[n/kpc]',&
					&'gradient Z direction[n/kpc]','gradient size[n/kpc]','average density at the sun[n]'
		WRITE(50,*)""

		DO e=min_cell_E,max_cell_E

			WRITE(*,*)min_cell_E,max_cell_E
			gradx=(density(x_max,y,z,e)-density(x_min,y,z,e))/(grid*2)

			grady=(density(x,y_max,z,e)-density(x,y_min,z,e))/(grid*2)

			gradz=(density(x,y,z_max,e)-density(x,y,z_min,e))/(grid*2)

			grad=SQRT(gradx**2+grady**2+gradz**2)
			WRITE(50,'(6f40.2)')10**(energy_grid*e),gradx,grady,gradz,grad,density(x,y,z,e)

		WRITE(50,*)""

		END DO

		CLOSE(50)
	END SUBROUTINE plot_anisotropy
END MODULE output
