doc/Calculate.cwg.R

In wendellopes/rvswf: RVSWF: Expansion in Vector Spherical Wave Functions written in R

#-------------------------------------------------------------------------------
# WAVE GUIDE PARAMETERS
#-------------------------------------------------------------------------------
rm(list=ls())
#-------------------------------------------------------------------------------
# Zeros of Bessel Functions
#-------------------------------------------------------------------------------
#Nth zero (lines) of the Mth Bessel functions (column)
ZJ <-matrix(c( 2.4048, 3.8317, 5.1356, 6.3802, 7.5883, 8.7715,
               5.5201, 7.0156, 8.4172, 9.7610,11.0647,12.3386,
               8.6537,10.1735,11.6198,13.0152,14.3725,15.7002,
              11.7915,13.3237,14.7960,16.2235,17.6160,18.9801,
              14.9309,16.4706,17.9598,19.4094,20.8269,22.2178),
     ncol=6,byrow=TRUE)
#-------------------------------------------------------------------------------
# Zeros of derivatives of Bessel Functions
#-------------------------------------------------------------------------------
#Nth zero (lines) of the derivative of Mth Bessel functions (column)
ZdJ<-matrix(c( 3.8317, 1.8412, 3.0542, 4.2012, 5.3175, 6.4156,
               7.0156, 5.3314, 6.7061, 8.0152, 9.2824,10.5199,
              10.1735, 8.5363, 9.9695,11.3459,12.6819,13.9872,
              13.3237,11.7060,13.1704,14.5858,15.9641,17.3128,
              16.4706,14.8636,16.3475,17.7887,19.1960,20.5755),
     ncol=6,byrow=TRUE)
#-------------------------------------------------------------------------------
# Wave Guide Parameters
#-------------------------------------------------------------------------------
# Basic Parameters
lambda=.5e-6            # Propagating wavelength
k=2*pi/lambda           # Propagating wavenumber
S=-1                    # Chirality of the beam
M<-5  # 0<=M<=5         # Order of the (derivative of the) Bessel function
N<-4  # 1<=N<=5         # Order of the zero of the derivative/Bessel function
TM<-TRUE                # Mode: TM (Bessel function), TE (derivative)
R<-5*lambda             # Radius of the waveguide
#-------------------------------------------------------------------------------
# MODES: M <- M+1 because J_0(x) ocupies column 1
# TE - ZdJ (derivative)
# TM - ZJ  (Bessel function)
#-------------------------------------------------------------------------------
if(TM){
   g<-ZJ[N,M+1]/R
   MODE<--1
   MODE.PWE<-4
}else{
   g<-ZdJ[N,M+1]/R
   MODE<-1
   MODE.PWE<-3
}
kz<-sqrt(k^2-g^2)
#-------------------------------------------------------------------------------
# Geometry of the calculations
#-------------------------------------------------------------------------------
NPX<-4*50+1
NPY<-4*60+1
NPZ<-4*20+1
dx<-2*R/(NPX-1)
dy<-2*R/(NPY-1)
dz<-2*R/(NPZ-1)
x<-seq(-R,R,by=dx)
y<-seq(-R,R,by=dy)
z<-seq(-R,R,by=dz)
#-------------------------------------------------------------------------------
lmax<- 2  # Number of partial waves to sum (for some l we have 2l+1 values of m)
#-------------------------------------------------------------------------------
# POSITION AT WHICH THE EXPANSION WILL BE PERFORMED  (REFERENCE SYSTEM)
#-------------------------------------------------------------------------------
# ARBITRARY
set.seed(512)
xo<-sample(x,1)
yo<-sample(y,1)
zo<-sample(z,1)
# FIXED
xo<-x[NPX%/%4+1]
yo<-y[NPY%/%4+1]
xo<-yo<-zo<-0
#-------------------------------------------------------------------------------
# CHANGE THE REFERENCE SYSTEM TO THE NEW POSITIONS
#-------------------------------------------------------------------------------
x<-x-(xo+dx/2)
y<-y-(yo+dy/2)
z<-z-(zo+dz/2)
z<-0;NPZ<-1
#-------------------------------------------------------------------------------
# BSC CALCULATIONS
#-------------------------------------------------------------------------------
CWG<-vwfd.cwg(TE=!TM,M,S,g,kz,x+xo,y+yo,z+zo)
BSC<-vswf.cwg(TM,g,kz,xo,yo,zo,lmax,m=M,s=S)
PWE<-vswf.pwe(k,x,y,z,lmax,BSC$GTE,BSC$GTM)
if(TM){ # TM implies Hz=0
   tez.CWG<-array(CWG$Ez,c(NPZ,NPY,NPX))[1,,]
   tez.PWE<-array(PWE$Ez,c(NPZ,NPY,NPX))[1,,]
}else{  # TE implies Ez=0
   thz.CWG<-array(CWG$Hz,c(NPZ,NPY,NPX))[1,,]
   thz.PWE<-array(PWE$Hz,c(NPZ,NPY,NPX))[1,,]
}
#-------------------------------------------------------------------------------
# IMAGE
#-------------------------------------------------------------------------------
source("plots.cwg.R")
#-------------------------------------------------------------------------------
#source("Image.R")
#u<-seq(0,2*pi,length.out=200)
#if(TM){
#   zl<-range(Re(tez.CWG))
#   x11()
#   Image((y+yo)/lambda,(x+xo)/lambda,z=Re(tez.CWG),nlevels=256,axes=TRUE,
#      color.palette=cm.colors,asp=1,#zlim=zl,
#      plot.axes={
#         axis(1);
#         axis(2);
#         abline(h=xo/lambda,v=yo/lambda,col='green')
#         lines(R/lambda*cbind(sin(u),cos(u)),col='red',lwd=2)
#      },
#      xlab=expression(y/lambda),ylab=expression(x/lambda))
#   x11()
#   Image((y+yo)/lambda,(x+xo)/lambda,z=Re(tez.PWE),nlevels=256,axes=TRUE,
#      color.palette=cm.colors,asp=1,#zlim=zl,
#      plot.axes={
#         axis(1);
#         axis(2);
#         abline(h=xo/lambda,v=yo/lambda,col='green')
#         lines(R/lambda*cbind(sin(u),cos(u)),col='red',lwd=2)
#      },
#      xlab=expression(y/lambda),ylab=expression(x/lambda))
#}else{
#   zl<-range(Re(thz.CWG))
#   x11()
#   Image((y+yo)/lambda,(x+xo)/lambda,z=Re(thz.CWG),nlevels=256,axes=TRUE,
#      color.palette=cm.colors,asp=1,#zlim=zl,
#      plot.axes={
#         axis(1);
#         axis(2);
#         abline(h=xo/lambda,v=yo/lambda,col='green')
#         lines(R/lambda*cbind(sin(u),cos(u)),col='red',lwd=2)
#         },
#      xlab=expression(y/lambda),ylab=expression(x/lambda))
#   x11()
#   Image((y+yo)/lambda,(x+xo)/lambda,z=Re(thz.PWE),nlevels=256,axes=TRUE,
#   color.palette=cm.colors,asp=1,#zlim=zl,
#      plot.axes={
#         axis(1);
#         axis(2);
#         abline(h=xo/lambda,v=yo/lambda,col='green')
#         lines(R/lambda*cbind(sin(u),cos(u)),col='red',lwd=2)
#      },
#      xlab=expression(y/lambda),ylab=expression(x/lambda))
#}
#