	program comparison
! Auhtor: Toni Pujol
! PLEASE OPEN def_sim_lpn.f90 FIRST, AND FOLLOW THE INSTRUCTIONS THERE.
! VERY IMPORTANT!! PLEASE CLOSE THAT PROJECT BEFORE COMPILING AND EXECUTING THIS PROGRAM (FILE-->CLOSE WORKSPLACE IF YOU USE FORTRAN POWERSTATION 4.0)
!
! Program to compare simulations with archaeological data
!
! This program computes the mean over all sites of the magnitude of /archaeological date - simulated date/ (table 1 in the paper) and shows it on the screen.
!
! This program also creates the output.dat file. The columns of the output.dat file are:
!		latitude (º), longitude (º), archaeological site date (cal yr BP), simulated site date (cal yr BP), distance site-origin (km)
! 		[origin=Jericho (for initial conditions 0 below) or Hemar (for initial conditions 1 or 2 below)]

	implicit none

	integer ::	index,indey,ix(20000),iy(20000),isimu,idumm
	integer ::  ipoin,jpoin,indad,indx,indy,indxm,indym
	integer ::  iskip
	integer ::  icondini

	character*10 :: fileinput
	
	real*8 :: ysim(200,200),varf,sddif,isitey,isitex
	real*8 :: x1(2000),x2(2000),y12(2000),xlat(2000),xlon(2000),dist(2000)
	real*8 :: ysim1,ysim2,ysim3,ysim4,x2ad,x2au,x1ad,x1au
	real*8 :: yinte,coeft,coefu,diffy,reldi
	real*8 :: varpi,xval,yval,coordx,coordy,coefc,coefn,coefr
	real*8 :: reflat1,reflat2,vorilat,vorilon,rho,theta
	real*8 :: vix,viy,datehemar,y12bp,sdiff
	real*8 :: xlatdeg(2000),xlondeg(2000),gennum,tgen

	write(*,*) 'Please enter 0, 1 or 2 for the initial conditions you simulated with sim_ind_micro.f90 (if you are not sure, this is parameter icondini in sim_ind_micro.inp):'
	write(*,*) ' '
	write(*,*) '0 = 87 CONTEXT sites, Jericho as origin (first row in Table 1 and Figs.2,3)'
	write(*,*) '1 = 16 PPNB/C sites, Hemar as origin (second row in Table 1)'
	write(*,*) '2 = 16 PPNB/C sites at 9000 cal yr BP (last row in Table 1 and Figs.7,8,B1,B2)'
	write(*,*) ' '

	read(*,*) icondini
	write(*,*) ' '

	if(icondini .eq. 0) then
		datehemar = 11863. ! date Jerico
	else if(icondini.eq. 1) then
		datehemar = 10418.	! date Hemar
	else if(icondini.eq. 2) then
		 datehemar = 9000.
	else
		write(*,*) 'write option 0, 1 or 2'
		stop
	end if

	write(*,*) 'Years per generation (e.g., 32)'

	read(*,*) tgen
	write(*,*) ' '

!
! Initialize at 0
!
	do ipoin = 1,20000

		ix(ipoin) = 0
		iy(ipoin) = 0

	end do

	do ipoin = 1,200

		do jpoin = 1,200

			ysim(ipoin,jpoin) = 0.

		end do

	end do

	do ipoin = 1,2000

		x1(ipoin) = 0.
		x2(ipoin) = 0.
		y12(ipoin) = 0.
		xlat(ipoin) = 0.
		xlon(ipoin) = 0.

	end do

!
! end initialization
!

!
! Input file (11) => simulations 
! .
!
!
	write(*,*) 'Name input file	(e.g. n150_0.dat)'
	read(*,*) fileinput

	isimu = 0

	open (11,file=fileinput,status='old')


		do while(.not.EOF(11))

		isimu = isimu + 1 
	
		read(11,*) vix,viy,varf

		ix(isimu) = int(vix)
		iy(isimu) = int(viy)

		if(varf.gt.240.) varf = -150.

		ysim(ix(isimu),iy(isimu)) = (varf - 1.)* tgen

	end do

	close (11)

	index = 0
	indey = 0

	do ipoin = 1, isimu

		if(ix(ipoin).gt.index) index = ix(ipoin)
		if(iy(ipoin).gt.indey) indey = iy(ipoin)

	end do

!
! Input files (12) => data FEPRE + CONTEXT
!

	varpi = 3.141592653590

	indad = 0

	if(icondini.eq.0) then
		open (12,file='lat_lon_bef_jerico_11863.dat',status='old')
	else if(icondini.eq.1) then
		open (12,file='lat_lon_calBP_distHemar.dat',status='old')
	else if(icondini.eq.2) then
		open (12,file='lat_lon_calBP_distHemar_all9000.dat',status='old')
	end if

	if(icondini.eq.2 .OR. icondini.eq.1) then
		do while (.not.EOF(12) .AND. indad.le.918) 

			indad = indad + 1

			read(12,*) xval,yval,y12bp,dist(indad)
			
			y12(indad) = datehemar - y12bp

			xlat(indad) = xval*varpi/180.
			xlon(indad) = yval*varpi/180.

			xlatdeg(indad) = xval
			xlondeg(indad) = yval

		end do
	else
		do while (.not.EOF(12) .AND. indad.le.989) 
	
			indad = indad + 1

			read(12,*) xval,yval,y12bp,dist(indad)
			
			y12(indad) = datehemar - y12bp

			xlat(indad) = xval*varpi/180.
			xlon(indad) = yval*varpi/180.

			xlatdeg(indad) = xval
			xlondeg(indad) = yval

		end do
	end if

	close (12)


!
! From lon - lat of data to X - Y on the grid
!
! Albers projection
!

! Reference latitudes 1 and 2

reflat1 = 30.*varpi/180.
reflat2 = 65.*varpi/180.

! latitude and longitude origin

vorilat = 31.17*varpi/180.
vorilon = 35.18*varpi/180.

! coefficients for the Albers projection

coefn = .5*(sin(reflat1) + sin(reflat2))
coefc = cos(reflat1)*cos(reflat1)+2.*coefn*sin(reflat1)
coefr = sqrt(coefc-2.*coefn*sin(vorilat))/coefn

! Projection on the plane

do ipoin = 1,indad

	theta = coefn*(xlon(ipoin) - vorilon)
	rho = sqrt(coefc - 2.*coefn*sin(xlat(ipoin)))/coefn

	coordx = rho*sin(theta)
	coordy = coefr - rho*cos(theta)

	x1(ipoin) = (coordx + 0.8632866)/.00784806 + 1.
	x2(ipoin) = (coordy + 0.16480926)/.00784806 + 1.

end do



!
! Output file (13)
!

	open (13,file='output.dat')

!
! Comparison
!

	sdiff = 0.
	sddif = 0.
	idumm = 0

	do ipoin = 1, indad

!
!  Find quadrant in the middle for x1
!		

		indx = 1

		do while (indx.lt.x1(ipoin) .AND. indx.lt.index)
		
			indx = indx + 1

		end do

		if(indx.lt.index) then
			indxm = indx - 1 
			x1ad = indxm
			x1au = indxm+1

			isitex = 1

		else 

!			write(*,*) indx,ipoin

			isitex = 0

		end if

!
!  Find quadrant in the middle for x1
!		
		indy = 1

		do while (indy.lt.x2(ipoin) .AND. indy.lt.indey)
		
			indy = indy + 1

		end do

		if(indy.lt.indey) then
			indym = indy-1
			x2ad = indym
			x2au = indym+1

			isitey = 1
		
		else

!			write(*,*) indy,ipoin,x2(ipoin),indey
	
			isitey = 0

		end if
!
! Bilineal interpolation (Numerical Recipes, pag. 117)
!
		
		if(isitex.eq.1 .AND. isitey.eq.1) then !!! go to boundaries

			ysim1 = ysim(indxm,indym)
			ysim2 =	ysim(indxm+1,indym)
			ysim3 = ysim(indxm+1,indym+1)
			ysim4 =	ysim(indxm,indym+1)

			iskip = 0

			if(ysim1.lt.-1. .and. ysim2.lt.-1. .and. ysim3.lt.-1. .and. ysim4.lt.-1.) then
				iskip = 1
			
			else if(ysim1.lt.-1. .and. ysim2.lt.-1. .and. ysim3.lt.-1. .and. ysim4.gt.-1.) then
				ysim1=ysim4
				ysim2=ysim4
				ysim3=ysim4
			
			else if(ysim1.lt.-1. .and. ysim2.lt.-1. .and. ysim3.gt.-1. .and. ysim4.gt.-1.) then
				ysim1=ysim4
				ysim2=ysim3

			else if(ysim1.lt.-1. .and. ysim2.gt.-1. .and. ysim3.lt.-1. .and. ysim4.gt.-1.) then
				ysim1=ysim4
				ysim3=ysim2

			else if(ysim1.gt.-1. .and. ysim2.lt.-1. .and. ysim3.lt.-1. .and. ysim4.gt.-1.) then
				ysim2=ysim1
				ysim3=ysim4

			else if(ysim1.lt.-1. .and. ysim2.lt.-1. .and. ysim3.gt.-1. .and. ysim4.lt.-1.) then
				ysim1=ysim3
				ysim2=ysim3
				ysim4=ysim3

			else if(ysim1.lt.-1. .and. ysim2.gt.-1. .and. ysim3.gt.-1. .and. ysim4.lt.-1.) then
				ysim1=ysim2
				ysim4=ysim3

			else if(ysim1.gt.-1. .and. ysim2.lt.-1. .and. ysim3.gt.-1. .and. ysim4.lt.-1.) then
				ysim4=ysim1
				ysim2=ysim3

			else if(ysim1.lt.-1. .and. ysim2.gt.-1. .and. ysim3.lt.-1. .and. ysim4.lt.-1.) then
				ysim1=ysim2
				ysim3=ysim2
				ysim4=ysim2

			else if(ysim1.gt.-1. .and. ysim2.gt.-1. .and. ysim3.lt.-1. .and. ysim4.lt.-1.) then
				ysim4=ysim1
				ysim3=ysim2

			else if(ysim1.gt.-1. .and. ysim2.lt.-1. .and. ysim3.lt.-1. .and. ysim4.lt.-1.) then
				ysim2=ysim1
				ysim3=ysim1
				ysim4=ysim1

			else if(ysim1.lt.-1. .and. ysim2.gt.-1. .and. ysim3.gt.-1. .and. ysim4.gt.-1.) then
				ysim1=ysim4

			else if(ysim1.gt.-1. .and. ysim2.lt.-1. .and. ysim3.gt.-1. .and. ysim4.gt.-1.) then
				ysim2=ysim3

			else if(ysim1.gt.-1. .and. ysim2.gt.-1. .and. ysim3.lt.-1. .and. ysim4.gt.-1.) then
				ysim3=ysim2

			else if(ysim1.gt.-1. .and. ysim2.gt.-1. .and. ysim3.gt.-1. .and. ysim4.lt.-1.) then
				ysim4=ysim1

			else if(ysim1.gt.-1. .and. ysim2.gt.-1. .and. ysim3.gt.-1. .and. ysim4.gt.-1.) then

			else

				write(*,*) ysim1,ysim2,ysim3,ysim4
				stop

			end if

			if(iskip.eq.0) then

				coeft = (x1(ipoin)-x1ad)/(x1au-x1ad)
				coefu = (x2(ipoin)-x2ad)/(x2au-x2ad)

				yinte = (1.-coeft)*(1.-coefu)*ysim1 + coeft*(1.-coefu)*ysim2 + coeft*coefu*ysim3 + (1.-coeft)*coefu*ysim4

! 
!
					idumm = idumm + 1

					diffy = abs(yinte-y12(ipoin))

					sdiff = sdiff + diffy

					sddif = sddif + diffy*diffy 

!					write(*,*) y12(ipoin)

					if(y12(ipoin).eq.0.) then
						reldi = 0.
					else
						reldi = diffy/(datehemar-y12(ipoin))

					end if

					if(ipoin.eq.908) yinte=0.	

					gennum = yinte/tgen

					write(13,400) xlatdeg(ipoin),xlondeg(ipoin),datehemar-y12(ipoin),datehemar-yinte,dist(ipoin)

			end if

		
		else if(isitex.eq.1. .AND. isitey.eq.0) then

			
			if(x2(ipoin).gt.indey-1 .AND. x2(ipoin).lt.indey) then
			
				indym = indey-1

				ysim1 = ysim(indxm,indym)
				ysim2 =	ysim(indxm+1,indym)
				ysim3 = ysim(indxm+1,indym+1)
				ysim4 =	ysim(indxm,indym+1)

				iskip = 0

				if(ysim1.lt.-1. .and. ysim2.lt.-1. .and. ysim3.lt.-1. .and. ysim4.lt.-1.) then

					iskip = 1

				else if(ysim1.lt.-1. .and. ysim2.lt.-1. .and. ysim3.lt.-1. .and. ysim4.gt.-1.) then
					ysim1=ysim4
					ysim2=ysim4
					ysim3=ysim4
				
				else if(ysim1.lt.-1. .and. ysim2.lt.-1. .and. ysim3.gt.-1. .and. ysim4.gt.-1.) then
					ysim1=ysim4
					ysim2=ysim3

				else if(ysim1.lt.-1. .and. ysim2.gt.-1. .and. ysim3.lt.-1. .and. ysim4.gt.-1.) then
					ysim1=ysim4
					ysim3=ysim2

				else if(ysim1.gt.-1. .and. ysim2.lt.-1. .and. ysim3.lt.-1. .and. ysim4.gt.-1.) then
					ysim2=ysim1
					ysim3=ysim4

				else if(ysim1.lt.-1. .and. ysim2.lt.-1. .and. ysim3.gt.-1. .and. ysim4.lt.-1.) then
					ysim1=ysim3
					ysim2=ysim3
					ysim4=ysim3

				else if(ysim1.lt.-1. .and. ysim2.gt.-1. .and. ysim3.gt.-1. .and. ysim4.lt.-1.) then
					ysim1=ysim2
					ysim4=ysim3

				else if(ysim1.gt.-1. .and. ysim2.lt.-1. .and. ysim3.gt.-1. .and. ysim4.lt.-1.) then
					ysim4=ysim1
					ysim2=ysim3

				else if(ysim1.lt.-1. .and. ysim2.gt.-1. .and. ysim3.lt.-1. .and. ysim4.lt.-1.) then
					ysim1=ysim2
					ysim3=ysim2
					ysim4=ysim2

				else if(ysim1.gt.-1. .and. ysim2.gt.-1. .and. ysim3.lt.-1. .and. ysim4.lt.-1.) then
					ysim4=ysim1
					ysim3=ysim2

				else if(ysim1.gt.-1. .and. ysim2.lt.-1. .and. ysim3.lt.-1. .and. ysim4.lt.-1.) then
					ysim2=ysim1
					ysim3=ysim1
					ysim4=ysim1

				else if(ysim1.lt.-1. .and. ysim2.gt.-1. .and. ysim3.gt.-1. .and. ysim4.gt.-1.) then
					ysim1=ysim4

				else if(ysim1.gt.-1. .and. ysim2.lt.-1. .and. ysim3.gt.-1. .and. ysim4.gt.-1.) then
					ysim2=ysim3

				else if(ysim1.gt.-1. .and. ysim2.gt.-1. .and. ysim3.lt.-1. .and. ysim4.gt.-1.) then
					ysim3=ysim2

				else if(ysim1.gt.-1. .and. ysim2.gt.-1. .and. ysim3.gt.-1. .and. ysim4.lt.-1.) then
					ysim4=ysim1

				else if(ysim1.gt.-1. .and. ysim2.gt.-1. .and. ysim3.gt.-1. .and. ysim4.gt.-1.) then

				else

					write(*,*) ysim1,ysim2,ysim3,ysim4
					stop

				end if

				if(iskip.eq.0) then

					coeft = (x1(ipoin)-x1ad)/(x1au-x1ad)
					coefu = (x2(ipoin)-x2ad)/(x2au-x2ad)

					yinte = (1.-coeft)*(1.-coefu)*ysim1 + coeft*(1.-coefu)*ysim2 + coeft*coefu*ysim3 + (1.-coeft)*coefu*ysim4

	!
						idumm = idumm + 1

						diffy = abs(yinte-y12(ipoin))

						sdiff = sdiff + diffy

						sddif = sddif + diffy*diffy 

						if(y12(ipoin).eq.0) then
							reldi = 0.
						else
							reldi = diffy/(datehemar-y12(ipoin))

						end if

						gennum = yinte/tgen

						write(13,400) xlatdeg(ipoin),xlondeg(ipoin),datehemar-y12(ipoin),datehemar-yinte,dist(ipoin)

				end if

			else if(x2(ipoin).gt.indey-1 .AND. x2(ipoin).lt.indey) then

				stop

			end if

		else if(isitex.eq.0. .AND. isitey.eq.1) then

			write(*,*) 'still not written'
			
			stop

		else if(isitex.eq.0. .AND. isitey.eq.0) then
		
			write(*,*) 'out of boundaries'		
		
			stop

		end if

	end do

	write(*,*) ' '
	write(*,*) 'Mean over all sites of the magnitude of /archaeological site date - simulation site date/ (table 1) = ',sdiff/idumm,' cal yr BP'
	write(*,*) ' '
	write(*,*) 'You can find the latitude, longitude, archaeological and simulated dates (cal yr BP), and km to Jericho (for initial conditions 0 above) or Hemar (for initial conditions 1 or 2) of all sites in file output.dat'
	write(*,*) ' '


	
400   format(5f15.7)
401   format(10f15.7)
402	  format(2f15.7)
	
	close(13)

	end
!
	