# Plot nonlinear study

library(RColorBrewer)

#load("track_64p_1024t_Q031_linA.RData")
load("tmp_octu_Q0618.RData")

Nk = length(T)
Np = T[[1]]$Np
Nt = T[[1]]$Nt
Rm = T[[1]]$Rm
RM = T[[1]]$RM
a = numeric(); for(i in 1:Nk) a = c(a, T[[i]]$A)
Am = min(a)
#AM = max(a)
AM = 16
rm(a)

matlab.colors =colorRampPalette(
  c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow",
    "#FF7F00", "red", "#7F0000"))
my.seq.palette = colorRampPalette( brewer.pal(9, "YlGnBu") )
#colo = rev(my.seq.palette(Np))
#colo = rev(topo.colors(Np))

wi = 8; he = 6; pt = 10; fo = "Times"; fn = "tmp_octu_Q0618"
sy = 16; sc = .2; LEGEND=TRUE; PHASE.SPACE=TRUE
SKIPt=1; SKIPp=1

dev.new(width=wi, height=he)
par(oma=rep(0, 4), mar=rep(0.1, 4), bg="black", xaxs="i", yaxs="i")

Nhor=4; Nver=3
layout(matrix(1:(Nhor*Nver), nrow=Nver, ncol=Nhor, byrow=TRUE))

for(i in 1:Nk){

  k = T[[i]]$k

  # sort by initial amplitude
  ia = sort.int(T[[i]]$A[, 1], index=TRUE)$ix
  ia = ia[seq(1, Np, by=SKIPp)]
  colo = matlab.colors(Np)
  colo = colo[seq(1, Np, by=SKIPp)]

  if(PHASE.SPACE){
#  x = T[[i]]$phi[ia, 1]
#  y = T[[i]]$A[ia, 1]
  x = T[[i]]$A[ia, 1] * cos(T[[i]]$phi[ia, 1])
  y = -T[[i]]$A[ia, 1] * sin(T[[i]]$phi[ia, 1])

  # plot initial values
  plot(x, y, axes=FALSE, xlab="", ylab="",
       xlim=c(-AM, AM), ylim=c(-AM, AM),
       pch=sy, cex=sc, col=colo)

  for(j in seq(2, Nt, by=SKIPt)){
#    x = T[[i]]$phi[ia, j]
#    y = T[[i]]$A[ia, j]
    x = T[[i]]$A[ia, j] * cos(T[[i]]$phi[ia, j])
    y = -T[[i]]$A[ia, j] * sin(T[[i]]$phi[ia, j])
    points(x, y, pch=sy, cex=sc, col=colo)
  }
  
  # show e-lens edges
#  symbols(x=c(0, 0), y=c(0, 0), circles=c(Rm, RM), inches=FALSE,
#          fg="white", add=TRUE)

} else {

  x = 1:Nt

  plot(c(1, Nt), c(Am, AM), axes=FALSE, xlab="", ylab="",
       type="n")

  ic = 0
  for(j in ia){
    y = T[[i]]$A[j, ]; ic = ic+1
    points(x, y, type="l", col=colo[ic])
  }
}

  if(LEGEND) legend("bottomright", bty="n",
                    legend=paste("k =", formatC(k, di=3)),
                    text.col="white", text.font=2, cex=1)
}

dev.copy2pdf(file=paste(fn, ".pdf", sep=""), pointsize=pt, family=fo)
dev2bitmap(file=paste(fn, ".png", sep=""), width=wi, height=he, res=300)