rm(list = ls())

##
## Create Figure 6 in manuscript
##


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


Sys.setenv(tz = "UTC")

all_vel_list = list()

dcme=2

vcme_vec = c(800,1000,1200)

list_time0 = list()

list_timef = list()

codeDir = getwd()

projectDir = dirname(codeDir)

plotsDir=paste0(projectDir,'/plots')

resultsDir=paste0(projectDir,'/results')

for (ivel in 1:3) {
	vcme = vcme_vec[ivel]

	resultsDir=paste0(projectDir,'/results')

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

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

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

	tracerSliceDir = paste0(masDir, "/dir_slice_tracers/")

	hdfDir = paste0(masDir, "/dir_2D/")

	sliceDir = paste0(masDir, "/dir_slice_1au/")

	setwd(projectDir)

	## MAS time units to sec
	
	MAS_to_km = 6.96 * 1e+05

	m_to_au = 149597870700

	one_au = m_to_au/(MAS_to_km)/1000

	MAS_to_km_sec = 481.3711

	mas_to_sec = 1445.87

	##
	## 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)
	cat("\n Processing Results From: ", masDir,'\n\n')
	#dumpDir = "dir_dump_times"

	setwd(timeDir)
	
	dump.files = list.files(path = ".", pattern = "dump_times")

	ndump = length(dump.files)

	list_hdf_dump_time = list_times = list_mas_df = list_date_hdf = list()

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

	}


	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)

		date_hdf = hdf_dump_time * mas_to_sec + time0

		# 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

	}


	## 
	## Cross Times
##


	# Now read the tracer files for each of the ndumps
	
	setwd(tracerDir)

	## Cross Time in days
	##

	cross_time = times[1:ndump]

	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 = rep(NA, nfiles)

		for (i in 1:nfiles) {
			fname = files[i]
			data = h5read(fname, "/Data")
			tr[i] = data[1, 1]
			H5close()
		}


		index = which(tr[] >= one_au)[1]

		cross_time[idir] = times[index]

		setwd("../")
	}


	dft = data.frame(dates = cross_time, variable = paste0("Realization", 1:ndump))

	pars = c("Vr", "Br", "Bt", "Bn", "B", "Temp", "np")

	##
	## Conversion factors from MAS Units
##

	## MAS velocity units to km/sec
	
	vel_to_km_sec = 481.3711

	## MAS Magnetic field units to nT
	
	b_to_nT = 2.2068908 * 1e-04 * 1e+09

	##  Electron density to 1/cm^3
	
	ne_to_den = 1e+08

	## Temperature to Kelvin
	
	Temp_to_K = 2.807067 * 1e+07

	##
	## Load MAS results
##

	setwd(sliceDir)
	cat("\n In sliceDir: ", sliceDir, '\n')
	list_mas_df = list()
	
	for (ifile in 1:ndump) {


		file = paste0("slice_1au_z", ifile, ".txt")

		mas_df = read.table(file = file, header = FALSE)

		# This is the order 
		# ${icount} ${val_br} ${val_bt} ${val_bp} ${val_vr} ${val_vt} ${val_vp} ${val_rho} ${val_t}

		colnames(mas_df) = c("r", "t", "p", "Br", "Bt", "Bp", "Vr", "Vt", "Vp", "np", "Temp")

		mas_df$r = mas_df$t = mas_df$p = NULL

		mas_df$B = sqrt(mas_df$Br^2 + mas_df$Bt^2 + mas_df$Bp^2) * b_to_nT

		mas_df$V = sqrt(mas_df$Vr^2 + mas_df$Vt^2 + mas_df$Vp^2) * vel_to_km_sec

		mas_df$Br = mas_df$Br * b_to_nT

		mas_df$Bt = mas_df$Bt * b_to_nT

		mas_df$Bp = mas_df$Bp * b_to_nT

		mas_df$Vr = mas_df$Vr * vel_to_km_sec

		mas_df$Vt = mas_df$Vt * vel_to_km_sec

		mas_df$Vp = mas_df$Vp * vel_to_km_sec

		mas_df$np = mas_df$np * ne_to_den

		mas_df$Temp = mas_df$Temp * Temp_to_K

		##
		## Add the UTC time to the MAS DF
##
		times = list_times[[ifile]]

		mas_df$ddate = times[1:length(mas_df$B)]

		list_mas_df[[ifile]] = mas_df

	}
	
	all_vel_list[[ivel]] = list_mas_df

	nfiles = 162
	
	times = list_times[[1]]

	list_time0[[ivel]] = round(times[1], 'hour')

	list_timef[[ivel]] = round(times[nfiles], 'hour')


}

## Save the list of data frames

file1 <- paste0('list_mas_df_vcme_',vcme,'_dcme_',dcme,'.rds')
saveRDS(file=file1, list_mas_df)

vars = c("V", "B", "Temp", "np")

labels = c("(a)", "(b)", "(c)", "(d)")

var.name = c("Speed [km/sec]", "Magnetic Field Magnitude [nT]", "Temperature [K]", "Density [1/cm^3]")


time_v500 = array(0, c(ndump, 3))

## Will need to interpolate the tracers


#times_1h = seq(from=time0, to=timef, by = '1 hour')

pl = list()

icount = 1

for (ivar in vars) {

	
	dfvel=list()
	ymin = 1e+06

	ymax = -1e+06
	ii=0
	for (ivel in 1:3) {
		DF = list()
		list_mas_df = all_vel_list[[ivel]]
		vcme = vcme_vec[ivel]
		times_1h = seq(from=list_time0[[ivel]], to=list_timef[[ivel]], by = '1 hour')
		for (ifile in 1:ndump) {
			ii = ii+1
			df = list_mas_df[[ifile]]

			y = df[1:nfiles, ivar]
			x = times[1:nfiles]
			yi = approx(x = as.numeric(x), y = y, xout = as.numeric(times_1h), rule = 2)$y

			DF[[ii]] = data.frame(dates = times_1h[1:nfiles], var = y)

			colnames(DF[[ii]])[2] <- paste0("VCME", ii)
			
			ymin = min(ymin, df[, ivar])

			ymax = max(ymax, df[, ivar])

			DF[[ii]] = reshape::melt(DF[[ii]], id.var='dates')
			if (ivar == 'V') {
				index = which(y > 500)[1]
				time_v500[ifile, ivel] = times_1h[index]
			}

		}
		dfvel[[ivel]] <- do.call("rbind", DF)
	}

	xlimits = c(floor_date(times_1h[1], unit = "hour"), ceiling_date(times_1h[length(times_1h)], unit = "hour"))

	

	ylimits = c(ymin, ymax)

	#colnames(DF[[ifile]])[2] <- (ndump + 1)


	pl[[icount]] = ggplot(data = dfvel[[1]], aes(x = dates, y = value, group = variable, alpha = 0.8)) + geom_line(color='green4') + geom_line(data = dfvel[[2]], aes(x = dates, y = value, group = variable, alpha=0.8), color = 'coral')+ geom_line(data = dfvel[[3]], aes(x = dates, y = value, group = variable, alpha=0.8), color = 'cornflowerblue')
	pl[[icount]] = pl[[icount]] + scale_x_datetime(name = "", date_breaks = "1 day", limits = xlimits, date_labels = "%b-%d-%y") 
	pl[[icount]] = pl[[icount]] + scale_y_continuous(name = var.name[icount], limits = ylimits) + theme(text = element_text(size = 12, color = "gray20", face = "italic"), axis.text.x = element_text(face = "plain", 
		angle = 60, size = 12, vjust = 0.4), axis.text.y = element_text(face = "plain", size = 12), legend.position = "none", legend.title = element_blank())
	pl[[icount]] = pl[[icount]] + ggtitle(var.name[icount]) + theme(plot.title = element_text(hjust = 0.5))

	pl[[icount]] = pl[[icount]] + geom_vline(xintercept = as.POSIXct(time0), color = "grey", linetype = "dashed")

	pl[[icount]] = pl[[icount]] + annotate('text', x = (xlimits[1]+3600), y = ymax, label = labels[icount], size = 6.0) 

	icount = icount + 1
}

margin = theme(plot.margin = unit(c(0.7,0.7,0.7,0.7), "cm"))

mp = grid.arrange(grobs = lapply(pl, "+", margin), nrow = 2, ncol = 2)

##


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

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

cat("\n\nSaving Figure 6 in:" ,filename,'\n\n')

# time_v500 is in seconds - remove the minimum values from it and convert to minutes

time_v500 = time_v500 - min(time_v500)
	
time_v500 = time_v500 / 60.0

colnames(time_v500) = c('v800', 'v1000', 'v1200')

for (ivel in 1:3) {
	cat('\nFor Velocity of: ', vcme_vec[ivel],'\n')
	cat("\n Range, Mean, Median, SD and 95 CI for Crossing Time: \n")
	cat(range(time_v500[,ivel]), mean(time_v500[,ivel]), median(time_v500[,ivel]), sd(time_v500[,ivel]), quantile(time_v500[,ivel], probs = c(0.05, 0.95)),'\n')
}

write.csv(time_v500, paste0(plotsDir,'/figure6_spread_stats.csv'))

setwd(codeDir)

##