rm(list = ls())

##
## Plot Br for 12 ADAPT maps at R = 1 Rs
##  Creates Figure 1
##

library(lubridate)
library(grid)
library(gridExtra)
library(zoo)
library(date)
library(rhdf5)
library(fields)
library(plot3D)


codeDir = getwd()

projectDir = dirname(codeDir)

resultsDir = paste0(projectDir,'/results/')

plotsDir=paste0(projectDir,'/plots')

colorTable<- designer.colors(51, c( "blue","white", "red") )

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

MAS_to_km_sec = 481.3711

MAS_to_nT = 2.2068908 * 1e-4 * 1e+9

MAS_to_Gauss = 2.2068908

##

var = 'br'
myr = 'r0'


dataDir = paste0(resultsDir, myr,"_files/")

setwd(dataDir)

saturate = FALSE
colorTable<- designer.colors(256, c( "blue","white", "red") )
if (var == 'br') {
		
	if (myr == 'r0') {
		zlim= c(-10.1,10.1)
		zmax = 10.0
		zmin =-10.0	
		saturate = TRUE	
	} else {
		zlim=c(-0.001, 0.001)
		zmax = zlim[2]
		zmin = zlim[1]
	}
	fac = MAS_to_Gauss
	colorTable<- designer.colors(256, c( "blue","white", "red") )
} else {
	zlim=c(250,660)
	zmax = 700.
	zmin=250.
	fac = MAS_to_km_sec
	colorTable<-designer.colors(256, c('darkblue','blue','green', 'yellow', 'red'))  # designer.colors(256, c("red", 'yellow', 'green', 'blue') )#designer.colors(256, rainbow(12))
}


files = list.files(path='.', pattern = paste0(var,'_',myr,'_'))

nfiles = length(files)

ninety = 90.0

DEGREES_TO_RADIANS = pi/180.0


filename = paste0(plotsDir,'/figure1.pdf')

pdf(file= filename, width = 12, height = 12)

if (var == 'br') {
	if (myr == 'r0') {
		par(mar = c(2.5, 2.5, 2, 2), mfrow = c(4, 3))
	} else {
		par(mar = c(2.5, 2.5, 2, 3), mfrow = c(4, 3))
	}
} else {
	par(mar = c(2.5, 2.5, 2, 2), mfrow = c(4, 3))
}



for (i in 1:nfiles) {

	fname = files[i]

	data = h5read(fname, "/Data")
	theta = h5read(fname, "/dim1")
	phi = h5read(fname, "/dim2")

	H5close()

	nt = length(theta)

	np = length(phi)

	# reverse so that north pole is at the top
	
	save_data = data

	data = data * 0

	for (j in 1:nt) {
		data[j, 1:np] = save_data[(nt - j + 1), 1:np]
	}

	theta = rev(theta)

	# move to degrees
	
	colat = theta/DEGREES_TO_RADIANS

	phi = phi/DEGREES_TO_RADIANS
	# move to lat -90 at the south pole and +90 at the north pole
	
	lat = ninety - colat

	data = data * fac

	if (saturate == TRUE) {
		for (j in 1:nt) {
			for (k in 1:np) {
				if (data[j, k] < 0) 
					data[j, k] = max(data[j, k], zmin)
				if (data[j, k] > 0) 
					data[j, k] = min(data[j, k], zmax)
			}
		}
	}

	tit = paste0("ADAPT REALIZATION # ", i)
	if( i%%3 == 0) {
	image2D(t(data),x=phi,y=lat, zlim = zlim, main = tit,xlab=expression(paste(Phi, " [deg]")),ylab=expression(paste(Lambda, " [deg]")), axes=F, xlim=range(phi), ylim=range(lat), col = colorTable)	
	} else {
	image2D(t(data),x=phi,y=lat, zlim = zlim, main = tit,xlab=expression(paste(Phi, " [deg]")),ylab=expression(paste(Lambda, " [deg]")), axes=F, xlim=range(phi), ylim=range(lat), col = colorTable, colkey = FALSE)		
	}

	axis(1, at = seq(from = 0, to = 360, by = 90), labels = seq(from = 0, to = 360, by = 90))
	axis(2, at = c(-90, -45, 0, 45, 90), labels = c(-90, -45, 0, 45, 90), las=2)

}

dev.off()

##
##