title "Waveguide"

select	
	
  errlim=1e-2
gridarc=5
coordinates
   cartesian3

variables

  Azr
  Azi
 ! Aur
 ! Aui
 ! Bur
 ! Bui

definitions


m=1e-2

D1=0.553*m
D=0.100*m
D0=0.117*m
!s=staged(8,8.25,8.5,8.75,9,9.25,9.5,9.75,10,10.25,10.5,10.75,11,11.25,11.5,11.75,12,12.25,12.5,12.75,13)
s=staged(9,9.1,9.2,9.3,9.4,9.5,9.6)
! s=9.35
x0=0.539*m
y0=m*2.283/2
b=1.01*m
h=0.554*m
a=2.283*m
l=40*m
lm=25*m
Eo=8.85e-12
Uo=4*pi*1e-7
C=1/sqrt(Uo*Eo)
f=s*1e9
w=2*pi*f

 G=0.8*ustep(x-lm)*(x-lm)
 K=1
 K1=2.9
 Km=k1
 Gm=0
 Cm=C/sqrt(1)
 DS= (pi/4)*(D0^2-D^2)
jx=0
jy=0
Kc=h*w/cm
jz=(1/(DS))*ustep(-(D/2)**2 + (x-x0)**2  + (y-y0)**2)*ustep((D0/2)**2 - (x-x0)**2  - (y-y0)**2)*ustep(h-z)*(sin((h-z)*w/cm)/sin((h)*w/cm))
!jz=(1/(DS))*(sin((h-z)*w/cm)/sin((h)*w/cm))
!jzi=0
!jzr=FIT(jz)
jzr=FIT(jz*sin(kc))
jzi=FIT(jz*cos(kc))






{ *************************************************************************************************************}

  Axr=0
  Axi=0
  Ayr=0
  Ayi=0

 Ar=vector(Axr,Ayr,Azr)
 Ai=vector(Axi,Ayi,Azi)

 Ui=-(1/(w*Uo*((Eo*k)^2+(G/w)^2)))*((G/w)*div(Ai)-Eo*k*div(Ar))
Ur=(1/(w*Uo*((Eo*k)^2+(G/w)^2)))*((G/w)*div(Ar)+Eo*k*div(Ai))
 Exr=-dx(Ur)+w*Axi
 Exi=-dx(Ui)-w*Axr
 Eyr=-dy(Ur)+w*Ayi
 Eyi=-dy(Ui)-w*Ayr
 Ezr=-dz(Ur)+w*Azi
 Ezi=-dz(Ui)-w*Azr
 Hxr=(1/Uo)*dy(Azr)
 Hxi=(1/Uo)*dy(Azi)
 Hyr=-(1/Uo)*dx(Azr)
 Hyi=-(1/Uo)*dx(Azi)

 DDPr= val(Ur,x0+D0/2+0.01*m,y0,0)-val(Ur,x0+(D1/2),y0,0)
 DDPi= val(Ui,x0+D0/2+0.01*m,y0,0)-val(Ui,x0+(D1/2),y0,0)

 Pr=-surf_integral(Ezr*Hyr+Ezi*Hyi,"F1")
 Pc=-surf_integral(Hyr*Ezi-Hyi*Ezr,"F1")

 Ir=-w*Azi*Jzr+w*Azr*Jzi
 Ii=w*Azi*Jzi-w*Azr*Jzr

 Ra=  Integral(Ir)
 Xa=  Integral(Ii)

Dn=(50+Ra)**2+Xa**2

Rho=(1/Dn)*sqrt((Ra**2-50**2+Xa*0**2)**2+4*(Xa*0*50)**2)

ROE=(1+Rho)/(1-Rho)

 !Ra= -Integral(Ezr*Jzr+Ezi*Jzi)
 !Xa=  Integral(Ezi*Jzr-Ezr*Jzi)
equations

     Azr: div(grad(Azr))+Uo*w^2*(Azr*Eo*k+(G/w)*Azi) =-Uo*Jzr
     Azi: div(grad(Azi))+Uo*w^2*(Azi*Eo*k-(G/w)*Azr) =-Uo*Jzi

 !    Azr: div(grad(Azr))+Uo*w^2*(Azr*Eo*k+(G/w)*Azi) =0
 !    Azi: div(grad(Azi))+Uo*w^2*(Azi*Eo*k-(G/w)*Azr) =0

 extrusion

 surface "Up" z=0
layer "metal"
 surface "med" z=h
layer "aire"
 surface "Down" z=b

boundaries
 region 1
!mesh_spacing = 0.005
 surface "up" natural(Azr)=0 natural(Azi)=0
 surface "down" natural(Azr)=0 natural(Azi)=0
  start(0,0)
  value(Azr)=0  value(Azi)=0
  line to (l,0)
  value(Azr)=0  value(Azi)=0
  line to (l,a)
  value(Azr)=0  value(Azi)=0
  line to (0,a)
  value(Azr)=0  value(Azi)=0
  line to finish

 feature "F1"
 start(7*m,0)
!mesh_spacing=0.5*m
 line to (7*m,a)


 region 2
  K=K1
! mesh_spacing=0.001
 surface "Up"  nobc(Azr) nobc(Azi)
 surface "med" nobc(Azr) nobc(Azi)
 surface "down" natural(Azr)=0 natural(Azi)=0
 start(x0-D1/2,y0)
 nobc(Azr) nobc(Azi)
 nobc(Azr) nobc(Azi)
 arc(center=x0,y0) angle=360

limited region 1 "F2"

 k=K1

 layer "metal" k=Km g=gm
 !surface "up" nobc(Azr) nobc(Azi)
 !surface "down" natural(Azr)=0 natural(Azi)=0
 !surface "med"  nobc(Azr) nobc(Azi)

 start(x0-D0/2,y0)
!mesh_spacing=0.0001
 layer "metal"   nobc(Azr) nobc(Azi)
 arc(center=x0,y0) angle=360



limited region 1 "F3"

 k=K1

  !layer "metal" void
 layer "metal" k=Km g=gm
 !surface "up" nobc(Azr) nobc(Azi)
 !surface "down" natural(Azr)=0 natural(Azi)=0
 !surface "med" nobc(Azr) nobc(Azi)

 start(x0-D/2,y0)
!mesh_spacing=0.00025
 !nobc(Azr) nobc(Azi)
 natural(Azr)=0 natural(Azi)=0
 arc(center=x0,y0) angle=360

 {

feature "Aq"
start(x0-D/2,y0)
mesh_spacing = 0.05*m
 arc(center=x0,y0) angle=360

 }



  monitors
  for cycle = 1
  contour(Azr^2+Azi^2)  on y=y0 zoom(0,0,2*Y0,h)
  contour(Azr^2+Azi^2)  on x=x0
!  contour(Ur)  on x=x0
!  contour(Ui)  on y=y0


  plots


 contour(Azr^2+Azi^2)  on x=x0
 vector(Eyr,Ezr)on x=x0
 contour(Azr^2+Azi^2) painted on y=y0
 contour(Azr^2+Azi^2) painted on z=h
 vector(Exr,Eyr) on z=0 zoom(x0-2*D0,y0-2*D0,4*D0,4*D0)
 elevation(Ezr*Hyr+Ezi*Hyi) from (5*m,y0,b/2) to (l,y0,b/2)
 contour(Ur^2+Ui^2) on z=0 zoom(x0-2*D0,y0-2*D0,4*D0,4*D0)
 grid(x,y) on z=0 zoom(x0-2*D0,y0-2*D0,4*D0,4*D0)


! transfer(Azr) file="Azr.dat"
! transfer(Azi) file="Azi.dat"



SUMMARY
REPORT  (DDPr^2 + DDPi^2)*Pr/(Pr^2 + Pc^2) as "R  "
REPORT  (DDPr^2 + DDPi^2)*Pc/(Pr^2 + Pc^2) as "X  "
!REPORT  val((jz^2)*(pi*(D/2)^2)^2,x0,y0,0.01*m)
REPORT  Ra as "Ra "
REPORT  Xa as "Xa "
REPORT  (Pr) as "Pr"
REPORT  (Pc) as "Pc"

HISTORIES

  HISTORY(Ra)
  HISTORY (Xa)
  HISTORY(Pr)
  HISTORY(Pc)
  HISTORY(ROE)
end