rm(list = ls())

##
## Create Figure 5 for Manuscript
##

library(lubridate)
library(forecast)
library(grid)
library(gridExtra)
library(zoo)
library(date)
library(ggplot2)
library(rhdf5)

Sys.setenv(tz = "UTC")


vcme=1000

dcme=2

codeDir = getwd()

projectDir = dirname(codeDir)

resultsDir=paste0(projectDir,'/results')

masDir = paste0(resultsDir,'/vcme',vcme,'_dcme',dcme)

timeDir = paste0(masDir, "/dir_dump_times/")

tracerDir = paste0(masDir, "/dir_tracers/")

plotsDir=paste0(projectDir,'/plots')


## MAS time units to sec

mas_to_sec = 1445.87 

MAS_to_km = 6.96 * 1e+05

m_to_au = 149597870700.

one_au = m_to_au / (MAS_to_km) / 1000.
##
## run id details
##

yyyy0 = 2019

mm0 = 2

dd0 = 13

hh0 = 0

mm20 = 0
 
time0  = ISOdatetime(year = yyyy0, month = mm0, day = dd0, hour = hh0, min = mm20, sec = 0, tz = "UTC")


##
## read the log file of the HDF dump times for this run
## 

##
## read the log file of the HDF dumps for this run
##

setwd(masDir)

subDir = "plots"

if(!file.exists(subDir)) {
        dir.create(subDir)
}

dumpDir='dir_dump_times'

setwd(dumpDir)

dump.files=list.files(path='.', pattern = 'dump_times')

ndump = length(dump.files)

list_hdf_dump_time = list_times = list_date_hdf = list_mas_df = list()

for (i in 1:ndump) {
	list_hdf_dump_time[[i]] = read.table(file=dump.files[i], header = FALSE)

}

# Process the times


for (i in 1:ndump) {
	##
	## need to substract relaxation time
	##
	hdf_dump_time = list_hdf_dump_time[[i]]
	
	hdf_dump_time = hdf_dump_time[, 1]
	
	deltat_ip = hdf_dump_time[1] #600.1360786600

	hdf_dump_time = hdf_dump_time - deltat_ip

	##
	## Convert to sec and add to time zero
	##

	date_hdf = hdf_dump_time * mas_to_sec + time0

	nt = length(hdf_dump_time)

	# convert to seconds
	
	hdf_dump_time = hdf_dump_time * mas_to_sec

	# These are the HDF times
	
	times = time0 + hdf_dump_time

	list_times[[i]] = times
	
	list_date_hdf[[i]] = date_hdf

}

rad2deg = 180.0/pi

deg2rad = 1.0/rad2deg

radCone = 5.0

radCone2 = radCone * radCone


# Now read the tracer files for each of the ndumps

setwd(tracerDir)

list_DF = list_cross = list_cross_cone = list()

## Cross Time in days
##

for (idir in 1:ndump) {

	dirname = paste0('tracers_wsa_z', idir)

	cat("Processing: ", dirname, '\n')
	
	setwd(dirname)

	times = list_date_hdf[[idir]]

	files = list.files(path = '.', pattern = 'tracers_pos')
	
	nfiles = length(files)
	
	fname = files[1]
	data = h5read(fname, "/Data")
	H5close()
	
	ntracers=dim(data)[1]
	
	tr = array(NA, c(nfiles, ntracers))
	
	if (idir == 1) {
		cross_time = rep(NA, ntracers)
	}
	
	for (i in 1:nfiles) {
		fname = files[i]
		data = h5read(fname, "/Data")
		tr[i, ] = data[,1]	
		H5close()
	}
	
	for (j in 1:ntracers) {
		
		index = which(tr[,j] >= one_au)[1]
		
		cross_time[j] = times[index] - times[1]
	}
	
	## Find the tracers that at time zero were within 5 degrees of the central particle.  
	
	if (idir == 1) {
		fname = files[1]
		data = h5read(fname, "/Data")

		t0 = data[1,2]
		p0 = data[1,3]
		
		ip_cone=rep(NA, ntracers)
		for (j in 1:ntracers) {
			rad2 = (t0-data[j,2])**2 + (p0-data[j,3])**2
			rad2 = rad2 * rad2deg * rad2deg
			
			if (rad2 <= radCone2) ip_cone[j] = j
		}
		ip_cone = ip_cone[!is.na(ip_cone)]
	}
	
	df = data.frame(dates = times[1:nfiles], tracers = tr)
	
	DF = reshape::melt(df, id.var = "dates")
	
	#DF$variable <- idir
	
	list_DF[[idir]] <- DF
	
	list_cross[[idir]] = cross_time
	
	list_cross_cone[[idir]] = cross_time[ip_cone]
	
	setwd("../")	
}

xmax = range(list_cross)[2]
xmax = round(xmax * 24.) + 2.
xmin = range(list_cross)[1]
xmin = round(xmin * 24 )-2.

ymax = 14
dy = 3

tracers_small_cone = tracers_large_cone = NULL
pl = list()

for (idir in 1:ndump) {
	
	xc = list_cross[[idir]] * 24.0
	
	xc_cone = list_cross_cone[[idir]] * 24.0
	
	tracers_small_cone = cbind(tracers_small_cone, as.numeric(unlist(xc_cone)))
	tracers_large_cone = cbind(tracers_large_cone, as.numeric(unlist(xc )))
	
	mean_cross_time = mean(xc)
	mean_corss_time = round(mean_cross_time, digits = 1)
	xc <- data.frame(xc)
	xc_cone <- data.frame(xc_cone)
	title = paste0("ADAPT MAP # ", idir)
	pl[[idir]] = ggplot() + geom_bar( data=xc, aes(xc, colour = 'cornflowerblue', alpha = 0.6),fill = 'cornflowerblue', alpha = 0.6)
	pl[[idir]] =pl[[idir]] + geom_bar(data=xc_cone, aes(xc_cone, colour = 'coral'),fill = 'coral', alpha = 0.8)
	pl[[idir]] = pl[[idir]] + theme(text = element_text(size = 12, color = "gray20", face = "italic"), axis.text.x = element_text(face = "plain", angle = 0, 
	 size = 12, vjust = 0.4), axis.text.y = element_text(face = "plain", size = 12),  legend.title = element_blank(), legend.position = "none")
	pl[[idir]] =  pl[[idir]]  + theme(plot.title = element_text(hjust = 0.5)) #+ ggtitle(title)
	pl[[idir]] = pl[[idir]] + scale_x_continuous(name = "Arrival Time [hours]", breaks= seq(from=24, to = 120, by = 12), limits = c(xmin, xmax))
	pl[[idir]] = pl[[idir]] + scale_y_continuous(name = "count", breaks= seq(from=0, to = ymax, by = dy), limits = c(0, ymax))
	pl[[idir]] = pl[[idir]] + annotate('text', x = (xmin+xmax)/2, y = ymax, label = title, size = 4.0) 
	
}

mp = grid.arrange(grobs = pl, nrow = 3, ncol = 4)

filename= paste0(plotsDir,"/figure5.pdf")

ggsave(file = filename, plot = mp, device = "pdf", width = 14, height = 10, units = "in")

cat("\n\n Saving Figure 5 in:" ,filename, '\n\n')

setwd(codeDir)

##
## Statistics on the tracers in the small and large cone areas
##


cat("\n\nFor Small (red bars) cone: \n\n")
cat("Min/Max Arrival Time: ", min(tracers_small_cone), max(tracers_small_cone), '\n\n')
cat("Mean Median and SD of Arrival Time: ", mean(tracers_small_cone), median(tracers_small_cone), sd(tracers_small_cone), '\n\n')
cat('Quantiles\n')
print(quantile(tracers_small_cone,probs=c(0,0.05,0.25,0.5,0.75,0.95,1.0)))
cat("\n\nFor Large (blue bars) cone: \n\n")
cat("Min/Max Arrival Time: ", min(tracers_large_cone), max(tracers_large_cone), '\n\n')
cat("Mean Median and SD of Arrival Time: ", mean(tracers_large_cone), median(tracers_large_cone), sd(tracers_large_cone), '\n\n')
cat('Quantiles\n')
print(quantile(tracers_large_cone,probs=c(0,0.05,0.25,0.5,0.75,0.95,1.0)))