rm(list = ls())

##
## Profiles of vr, |B|, T and number density for 12 different realizations at 1AU
## Figure 3 in manuscript
##

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

Sys.setenv(tz = "UTC")

vcme=1200

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/")

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

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

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

plotsDir=paste0(projectDir,'/plots')

## 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)

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_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 = 1.e+08

## Temperature to Kelvin

Temp_to_K = 2.807067 * 1.e+07

##
## Load MAS results
##

setwd(sliceDir) 

for (ifile in 1:ndump) {


	file = paste0("slice_1au_z", ifile, ".txt")
    cat("Processing File:", file,"\n")
	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
	
}

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]")


## Will need to interpolate the tracers

nfiles = 162

times = list_times[[1]]

time0 = round(times[1], 'hour')

timef = round(times[nfiles], 'hour')

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

pl = list()

icount = 1

for (ivar in vars) {

	DF = list()

	ymin = 1e+06

	ymax = -1e+06
	cat("Processing Var", ivar, "\n")
	for (ifile in 1:ndump) {

		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[[ifile]] = data.frame(dates = times_1h[1:nfiles], var = y)

		colnames(DF[[ifile]])[2] <- paste0("Realization",ifile)

		ymin = min(ymin, df[, ivar])

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

		DF[[ifile]] = reshape::melt(DF[[ifile]], id.var = "dates")

		imax = which.max(y)

	}

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

	df <- do.call("rbind", DF)

	ylimits = c(ymin, ymax)
	
	colnames(DF[[ifile]])[2] <- (ndump + 1)

	
	pl[[icount]] = ggplot(data = df, aes(x = dates, y = value, group = variable, color = variable)) + geom_line()
	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))
	#if (ivar == "V" || ivar == "B") 
		#pl[[icount]] = pl[[icount]] + geom_line(data = dfc, aes(x = dates, y = value, group = variable, colour = value), size = 2)
	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) 
	
	
	#pl[[icount]] = pl[[icount]] + geom_vline(data = dft, aes(xintercept=dates, group = variable, color = variable), linetype='dashed')
	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,"/figure3.pdf")

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

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

setwd(codeDir)