CG  Version 25
annulusEigenFunction.h
Go to the documentation of this file.
1 //==================================================================================================
2 // Evaluate the annulus eigenfunction or it's error
3 //
4 // OPTION: solution, error
5 //==================================================================================================
6 #beginMacro annulusEigenFunction(OPTION)
7 
8 // printF(" I1.getBase(),uLocal.getBase(0),I1.getBound(),uLocal.getBound(0)=%i %i %i %i \n",
9 // I1.getBase(),uLocal.getBase(0),I1.getBound(),uLocal.getBound(0));
10 
11  Index J1 = Range(max(I1.getBase(),uLocal.getBase(0)),min(I1.getBound(),uLocal.getBound(0)));
12  Index J2 = Range(max(I2.getBase(),uLocal.getBase(1)),min(I2.getBound(),uLocal.getBound(1)));
13  Index J3 = Range(max(I3.getBase(),uLocal.getBase(2)),min(I3.getBound(),uLocal.getBound(2)));
14 
15  if( numberOfDimensions==2 )
16  {
17 #include "besselPrimeZeros.h"
18 
19  const int n = int(initialConditionParameters[0]+.5); // angular number, n=0,1,... --> Jn(omega*r)
20  const int m = int(initialConditionParameters[1]+.5); // radial number m=0,...
21 
22  assert( m<mdbpz && n<ndbpz );
23 
24  real omega = besselPrimeZeros[n][m]; // m'th zero of Jn' (excluding r=0 for J0)
25 
26 // printF("Annulus: Bessel function solution: n=%i, m=%i, omega=%e (c=%8.2e)\n",n,m,omega,c);
27 
28  const real eps=sqrt(REAL_EPSILON);
29 
30  real np1Factorial=1.;
31  for( int k=2; k<=n+1; k++ )
32  np1Factorial*=k; // (n+1)!
33 
34  int i1,i2,i3;
35  real r,gr,xd,yd,zd,bj,bjp,rx,ry,theta,thetax,thetay;
36  real cosTheta,sinTheta,bjThetax,bjThetay,uex,uey,cosn,sinn;
37 
38 
39  FOR_3D(i1,i2,i3,J1,J2,J3)
40  {
41  xd=X(i1,i2,i3,0);
42  yd=X(i1,i2,i3,1);
43  #If #DIM == "3"
44  zd=X(i1,i2,i3,2)
45  sinPhi=sin(Pi*zd)
46  cosPhi=cos(Pi*zd)
47  #End
48  r = sqrt(xd*xd+yd*yd);
49  theta=atan2(yd,xd);
50  // if( theta<0. ) theta+=2.*Pi;
51  cosTheta=cos(theta);
52  sinTheta=sin(theta);
53 
54  cosn=cos(n*theta);
55  sinn=sin(n*theta);
56 
57  gr=omega*r;
58 
59  rx = cosTheta; // x/r
60  ry = sinTheta; // y/r
61 
62  bj=jn(n,gr); // Bessel function J of order n
63 
64 
65  if( gr>eps ) // need asymptotic expansion for small gr ??
66  {
67  bjp = -jn(n+1,gr) + n*bj/gr; // from the recursion relation for Jn'
68  thetay= cosTheta/r;
69  thetax=-sinTheta/r;
70 
71  uex = (1./omega)*(omega*ry*bjp*cosn -n*bj*thetay*sinn);
72  uey = -(1./omega)*(omega*rx*bjp*cosn -n*bj*thetax*sinn);
73 
74 
75  }
76  else
77  {
78  // Jn(z) = (.5*z)^n *( 1 - (z*z/4)/(n+1)! + ..
79 
80 
81  // At r=0 all the Jn'(0) are zero except for n=1
82  // bjp = n==1 ? 1./2. : 0.;
83  bjp = n==0 ? 0. : pow(.5,double(n))*pow(gr,n-1.)*( 1. - (gr*gr)/(4.*np1Factorial) );
84 
85  // bj/r = omega*bjp at r=0
86  bjThetay= omega*bjp*cosTheta;
87  bjThetax=-omega*bjp*sinTheta;
88 
89  uex = (1./omega)*(omega*ry*bjp*cosn -n*bjThetay*sinn); // Ex.t = Hz.y
90  uey = -(1./omega)*(omega*rx*bjp*cosn -n*bjThetax*sinn); // Ey.t = - Hz.x
91 
92  }
93 
94  real sint = sin(omega*t), cost = cos(omega*t);
95 
96  #If #OPTION == "solution"
97  UHZ(i1,i2,i3) = bj*cosn*cost;
98  UEX(i1,i2,i3) = uex*sint; // Ex.t = Hz.y
99  UEY(i1,i2,i3) = uey*sint; // Ey.t = - Hz.x
100 
101  if( method==nfdtd )
102  {
103  UMHZ(i1,i2,i3) = bj*cosn*cos(omega*(t-dt));
104  UMEX(i1,i2,i3) = uex*sin(omega*(t-dt));
105  UMEY(i1,i2,i3) = uey*sin(omega*(t-dt));
106  }
107  else if( method==sosup )
108  {
109  uLocal(i1,i2,i3,hzt) = -omega*bj*cosn*sint;
110  uLocal(i1,i2,i3,ext) = omega*uex*cost;
111  uLocal(i1,i2,i3,eyt) = omega*uey*cost;
112  }
113 
114 
115  #Elif #OPTION == "error"
116  ERRHZ(i1,i2,i3) = UHZ(i1,i2,i3) -bj*cosn*cos(omega*t);
117  ERREX(i1,i2,i3) = UEX(i1,i2,i3) - uex*sin(omega*t); // Ex.t = Hz.y
118  ERREY(i1,i2,i3) = UEY(i1,i2,i3) - uey*sin(omega*t); // Ey.t = - Hz.x
119  if( method==sosup )
120  {
121  errLocal(i1,i2,i3,hzt) = uLocal(i1,i2,i3,hzt) + omega*bj*cosn*sint;
122  errLocal(i1,i2,i3,ext) = uLocal(i1,i2,i3,ext) - omega*uex*cost;
123  errLocal(i1,i2,i3,eyt) = uLocal(i1,i2,i3,eyt) - omega*uey*cost;
124  }
125  #Else
126  Overture::abort("error");
127  #End
128 
129  }
130  }
131  else /* 3D */
132  {
133 #include "besselZeros.h"
134 
135  const real cylinderLength=cylinderAxisEnd-cylinderAxisStart;
136 
137  const int n = int(initialConditionParameters[0]+.5); // angular number, n=0,1,... --> Jn(omega*r)
138  const int m = int(initialConditionParameters[1]+.5); // radial number m=0,...
139  const int k = int(initialConditionParameters[2]+.5); // axial number k=1,2,3,...
140 
141  assert( m<mdbz && n<ndbz );
142 
143  real lambda = besselZeros[n][m]; // m'th zero of Jn (excluding r=0 for J0)
144  real omega = sqrt( SQR(k*Pi/cylinderLength) + lambda*lambda );
145 
146  printF("***Cylinder: Bessel function soln: n=%i, m=%i, k=%i, lambda=%e, omega=%e (c=%8.2e) [za,zb]=[%4.2f,%4.2f]\n",
147  n,m,k,lambda,omega,c,cylinderAxisStart,cylinderAxisEnd);
148 
149  const real eps=sqrt(REAL_EPSILON);
150 
151  real np1Factorial=1.;
152  for( int k=2; k<=n+1; k++ )
153  np1Factorial*=k; // (n+1)!
154 
155  int i1,i2,i3;
156  real r,gr,xd,yd,zd,bj,bjp,rx,ry,theta,thetax,thetay;
157  real cosTheta,sinTheta,bjThetax,bjThetay,uex,uey,cosn,sinn,sinkz,coskz,cost,sint;
158 
159  FOR_3D(i1,i2,i3,J1,J2,J3)
160  {
161  xd=X(i1,i2,i3,0);
162  yd=X(i1,i2,i3,1);
163  zd=(X(i1,i2,i3,2)-cylinderAxisStart)/cylinderLength; // *wdh* 040626 -- allow for any length
164 
165  sinkz=sin(Pi*k*zd);
166  coskz=cos(Pi*k*zd);
167 
168  r = sqrt(xd*xd+yd*yd);
169  theta=atan2(yd,xd);
170 
171  cosTheta=cos(theta);
172  sinTheta=sin(theta);
173 
174  cosn=cos(n*theta);
175  sinn=sin(n*theta);
176 
177  cost=cos(omega*t);
178 
179  gr=lambda*r;
180 
181  rx = cosTheta; // x/r
182  ry = sinTheta; // y/r
183 
184  bj=jn(n,gr); // Bessel function J of order n
185 
186 
187  if( gr>eps ) // need asymptotic expansion for small gr ??
188  {
189  bjp = -jn(n+1,gr) + n*bj/gr; // from the recursion relation for Jn'
190  thetay= cosTheta/r;
191  thetax=-sinTheta/r;
192 
193  uex = -(k*Pi/(cylinderLength*lambda*lambda))*( lambda*rx*bjp*cosn - n*bj*thetax*sinn );
194  uey = -(k*Pi/(cylinderLength*lambda*lambda))*( lambda*ry*bjp*cosn - n*bj*thetay*sinn );
195 
196 
197  }
198  else
199  {
200  // Jn(z) = (.5*z)^n *( 1 - (z*z/4)/(n+1)! + ..
201 
202 
203  // At r=0 all the Jn'(0) are zero except for n=1
204  // bjp = n==1 ? 1./2. : 0.;
205  bjp = n==0 ? 0. : pow(.5,double(n))*pow(gr,n-1.)*( 1. - (gr*gr)/(4.*np1Factorial) );
206 
207  // bj/r = lambda*bjp at r=0
208  bjThetay= lambda*bjp*cosTheta;
209  bjThetax=-lambda*bjp*sinTheta;
210 
211  uex = -(k*Pi/(cylinderLength*lambda*lambda))*( lambda*rx*bjp*cosn -n*bjThetax*sinn); // Ex.t = Hz.y
212  uey = -(k*Pi/(cylinderLength*lambda*lambda))*( lambda*ry*bjp*cosn -n*bjThetay*sinn); // Ey.t = - Hz.x
213 
214  }
215 
216  #If #OPTION == "solution"
217 
218  UEX(i1,i2,i3) = uex*sinkz*cost;
219  UEY(i1,i2,i3) = uey*sinkz*cost;
220  UEZ(i1,i2,i3) = bj*cosn*coskz*cost;
221 
222  cost=cos(omega*(t-dt));
223  UMEX(i1,i2,i3) = uex*sinkz*cost;
224  UMEY(i1,i2,i3) = uey*sinkz*cost;
225  UMEZ(i1,i2,i3) = bj*cosn*coskz*cost;
226 
227 
228  #Elif #OPTION == "error"
229 
230 
231  ERREX(i1,i2,i3) = UEX(i1,i2,i3) - uex*sinkz*cost;
232  ERREY(i1,i2,i3) = UEY(i1,i2,i3) - uey*sinkz*cost;
233  ERREZ(i1,i2,i3) = UEZ(i1,i2,i3) - bj*cosn*coskz*cost;
234 
235  #Else
236  Overture::abort("error");
237  #End
238 
239  }
240  }
241 
242 #endMacro
Generated on Fri Jan 4 2013 10:19:08 for CG by   doxygen 1.8.3