## ----Loading_library,include=FALSE--------------------------------------------
knitr::opts_chunk$set(echo = TRUE, fig.align = 'center')
rm(list=ls())
library(RGeostats)
library(rgl)
library(knitr)
knit_hooks$set(webgl = hook_webgl)
variety.define(flag.sphere=TRUE)
constant.define("nproc",NA)
flag.plot = FALSE


## ----Creating_Grid------------------------------------------------------------
grid = db.create(nx=c(361,181),x0=c(0,90),dx=c(1,1))


## ----Creating_Model-----------------------------------------------------------
model = model.create(vartype="Cubic",range=0.5)


## ----Defining_Parameters------------------------------------------------------
nbfdef = 10000


## ----Spectrum_for_Chentsov_Covariance-----------------------------------------
set.keypair("Simsph_Shunt",1)
dummy = simsph(grid,model,special=1,verbose=TRUE)
freqs = get.keypair("Simsph_Spectrum_Frequencies")
plot(freqs, main="Chentsov Covariance Spectrum",
     type="l",ylab="Frequency")
plot(freqs, type="b", xlim=c(1,20), 
     main="Chentsov Covariance Spectrum (zoom)",ylab="Frequency")


## ----Spectrum_for_Exponential_Covariance--------------------------------------
set.keypair("Simsph_Shunt",1)
dummy = simsph(grid,model,special=2,verbose=TRUE)
freqs = get.keypair("Simsph_Spectrum_Frequencies")
plot(freqs, main="Exponential Covariance Spectrum",
     type="l",ylab="Frequency")
plot(freqs, type="b", xlim=c(1,20), 
     main="Exponential Covariance Spectrum (zoom)",ylab="Frequency")


## ----Spectrum_for_Standard_Model----------------------------------------------
set.keypair("Simsph_Shunt",1)
dummy = simsph(grid,model,special=0,verbose=TRUE)
freqs = get.keypair("Simsph_Spectrum_Frequencies")
plot(freqs, main="Standard Model Spectrum",type="l",ylab="Frequency")
plot(freqs, type="b", xlim=c(1,50), 
     main="Standard Model Spectrum (zoom)",ylab="Frequency")


## ----Function_ to_display_Spectrum_in_degrees---------------------------------
display.spectrum.degrees <- function(special=0,model,nbf=nbfdef,
                                     dslice=20,title="Undefined")
{
  set.keypair("Simsph_Shunt",2)
  dummy     = simsph(grid,model,special=special,verbose=FALSE,nbf=nbf)
  degrees   = get.keypair("Simsph_Array_Degree")
  freqs     = get.keypair("Simsph_Spectrum_Frequencies")
  fmax      = range(freqs)[2]
  dmax      = range(degrees)[2]
  h.result  = hist(degrees,breaks=seq(-0.5,0.5+dmax),plot=FALSE)
  nval      = length(h.result$mids)
  plot(h.result$mids,h.result$density,xlim=c(0,dslice),ylim=c(0,fmax),
       type="b",col="black",pch=19,
       xlab="Degree N",ylab="Frequency",main=title)
  lines(x=h.result$mids,y=freqs[1:nval],add=T,
        type="b",col="red",pch=19)
  set.keypair("Simsph_Shunt",0)
  invisible()
}


## ----Distribution_of_degrees_for_Chentsov_Covariance--------------------------
display.spectrum.degrees(special=1,model,
                         title="Chentsov Covariance Spectrum")


## ----Distribution_of_degrees_for_Exponential_Covariance-----------------------
display.spectrum.degrees(special=2,model,
                         title="Exponential Covariance Spectrum")


## ----Distribution_of_degrees_for_Standard_Model-------------------------------
display.spectrum.degrees(special=0,model,dslice=40,
                         title="Standard Model Spectrum")


## ----Function_to_display_Spectrum---------------------------------------------
display.spectrum.order <- function(special=0,model,nbf=nbfdef,
                                   title="Undefined")
{
  set.keypair("Simsph_Shunt",2)
  dummy     = simsph(grid,model,special=special,verbose=TRUE,nbf=nbf)
  degrees   = get.keypair("Simsph_Array_Degree")
  orders    = get.keypair("Simsph_Array_Order")
  freqs     = get.keypair("Simsph_Spectrum_Frequencies")
  dmax      = range(degrees)[2]
  fmax      = range(freqs)[2]
  h.result  = hist(orders,breaks=seq(-0.5-dmax,0.5+dmax),plot=FALSE)
  nval      = floor(length(h.result$mids) / 2) + 1
  plot(h.result$mids,h.result$density,xlim=c(-20,20),
       type="b",col="black",pch=19,
       xlab="Order K",ylab="Frequency",main=title)
  set.keypair("Simsph_Shunt",0)
  invisible()
}


## ----Distribution_of_orders_for_Chentsov_Covariance---------------------------
display.spectrum.order(special=1,model,
                      title="Chentsov Covariance Spectrum")


## ----Distribution_of_orders_for_Exponential_Covariance------------------------
display.spectrum.order(special=2,model,
                      title="Exponential Covariance Spectrum")


## ----Distribution_of_orders_for_Standard_Model--------------------------------
display.spectrum.order(special=0,model,
                      title="Standard Model Spectrum")


## ----Function_to_display_Spectrum_Phase---------------------------------------
display.spectrum.phase <- function(special=0,model,nbf=nbfdef,
                                   title="Undefined")
{
  set.keypair("Simsph_Shunt",2)
  dummy     = simsph(grid,model,special=special,verbose=FALSE,nbf=nbf)
  phases    = get.keypair("Simsph_Array_Phase")
  h.result  = hist(phases,breaks=100,
                   main="Histogram of phases",xlab="")
  set.keypair("Simsph_Shunt",0)
  invisible()
}


## ----Distribution_of_phases_for_Chentsov_Covariance---------------------------
display.spectrum.phase(special=1,model,
                      title="Chentsov Covariance Spectrum")


## ----Testing_Legendre_Function------------------------------------------------
set.keypair("Simsph_Test",1)
set.keypair("Simsph_Test_Degree",1)
set.keypair("Simsph_Test_Order",0)
set.keypair("Simsph_Test_Phase",0)
a = simsph(grid,model,special=1)
if (flag.plot) plot(a)


## ----Simulation_on_Sphere_with_Cubic_Model------------------------------------
model = model.create(vartype="Cubic",range=0.5,sill=1)
a = simsph(grid,model,nbf=1000)
if (flag.plot) plot(a)


## ----Simulation_on_Sphere_with_Matern_Model-----------------------------------
model = model.create(vartype="J-Bessel",range=0.03,param=0.01)
a = simsph(grid,model,nbf=1000)
if (flag.plot) plot(a)