Oblic emission and rotation troubles

How to use openEMS. Discussion on examples, tutorials etc

Moderator: thorsten

Youry
Posts: 51
Joined: Wed 22 Feb 2012, 12:50

Oblic emission and rotation troubles

Post by Youry » Mon 26 Mar 2012, 07:50

Dears,

I am trying to generate an oblic waves... I mean like a propagation as (1 1 0) direction.

Is their is an easy way to do this ? (because the permittivity geometry will become difficult to define if not...).

To do this I simply first imagine a rotation of the emission box...
Like :

Code: Select all

CSX = AddBox(CSX,'excitation',0,[10 15 10],[20 25 10],'Transform', {Rotate_Z, pi/4});
But this crash (not crash simply return error message)

Then I try more basic style, like rotate a metal cube... Also return error...

Code: Select all

% define the excitation
CSX = AddExcitation(CSX,'excitation',0,[0 1 0]); 
CSX = AddBox(CSX,'excitation',0,[10 15 10],[20 25 10]);

% geometry
CSX = AddMetal(CSX,'metal');
CSX = AddBox(CSX,'metal',10,[20 20 20],[22 22 22],'Transform', {Rotate_Z, pi/4});   [b]%This line bug[/b]
%CSX = AddBox(CSX,'metal',10,[20 20 20],[22 22 22]);                                             [b]%This line work[/b]


%Probing
CSX = AddDump(CSX,'Probe','DumpMode',0);
CSX = AddBox(CSX, 'Probe', 0, [0 4 0],[30 4 30]);
Did I forget somethings, because it is very simple thing, so I don't really understand... I also try :
AddBox(CSX,'metal',10,[20 20 20],[22 22 22],'Transform', {[0 0 1], pi/4});
AddBox(CSX,'metal',10,[20 20 20],[22 22 22],'Transform', {3, pi/4});

Do you see my 'Transform' problem ?
Is it the good way to make oblic emission ? If not... Is it possible to change the sinus phase, depending on the position in the box emission area ?

Thanks for advices, it will very very simplify my permittivity definition...

thorsten
Posts: 1438
Joined: Mon 27 Jun 2011, 12:26

Re: Oblic emission and rotation troubles

Post by thorsten » Mon 26 Mar 2012, 08:57

Hi,

what kind of an error do you get for the transformed boxes? I don't see any error in your first box definition. Did AppCSXCAD show it correctly (in 3D view mode)?
But still I'm not sure that this would do what you want anyway. But I must confess, I don't understand you goal.

Wouldn't it be easier to define the excitation in a plane and transform everything else? It is much easier to define a geometry dependent permittivity instead of such a strange excitation...

regards
Thorsten

Youry
Posts: 51
Joined: Wed 22 Feb 2012, 12:50

Re: Oblic emission and rotation troubles

Post by Youry » Mon 26 Mar 2012, 10:46

This one works :

Code: Select all

CSX = AddExcitation(CSX,'excitation',0,[0 1 0]); 
CSX = AddBox(CSX,'excitation',0,[0 0 10],[30 30 10]);
This one bugs :

Code: Select all

CSX = AddExcitation(CSX,'excitation',0,[0 1 0]); 
CSX = AddBox(CSX,'excitation',0,[0 0 10],[30 30 10],'Transform', {Rotate_Z, pi/4});
The error is simply that Rotate_Z is undefined.
There is also error evaluating list element number 1 and 7 (but it come from the unknown Rotate_Z I think).
That with I supposed that Rotate_Z was a vector and I try [0 0 1] in the 'Transform' argument but it seems it is not that...
AppCSXCAD is not called, the programme failed before it asked me if I want to delete all the previous result...

About the "strange" emission. I always try to define my permitivity geometry as closer as the geometry.
For example, if I want to simulate a simple 1D grating (rectangular shape), with incident waves at 30 degrees. If I decide to move the grating, all the straight line shape of the grating becomes a kind of stair, and I don't like that. It is probably not a problem if the grid-cell is small enough, but my computer is very limited. I feel that putting oblic wave excitation was good idea, it becomes very easy to change angle -instead of all the geometry- and my rectangular shape can fit exactly to the cartesian grid...
Maybe my mistake, I don't have FDTD experience in simulation result and precision, so I don't know how much this stairs (instead of line) affect the result.

Anyways, rotation is style not good, and I don't know why.

Moreover, is it possible to put all the permittivity geometry in 1 object, define a center, and rotate it all ?
(for example I have several long rectangular box (the grating), and rotate all in 1 time, not sure if it is clear ?)
Like this I can do what you propose, easily, first define the permittivity in the appropriate position with the grid, and then arbitrary rotated it... And keep propagation along x, y or z...
Last edited by Youry on Mon 26 Mar 2012, 11:03, edited 1 time in total.

thorsten
Posts: 1438
Joined: Mon 27 Jun 2011, 12:26

Re: Oblic emission and rotation troubles

Post by thorsten » Mon 26 Mar 2012, 16:51

Hi,

sorry I missed it, but Rotate_Z is a string so you need the ' '

Code: Select all

CSX = AddBox(CSX,'excitation',0,[0 0 10],[30 30 10],'Transform', {Rotate_Z, pi/4});  %wrong
CSX = AddBox(CSX,'excitation',0,[0 0 10],[30 30 10],'Transform', {'Rotate_Z', pi/4}); %correct
Every object needs its own transformation, but just define a matlab variable and do copy and paste...
Then you only need to change alpha once or even do a sweep of alpha
e.g.

Code: Select all

alpha = pi/4;
CSX = AddBox(CSX,'excitation',0,[0 0 10],[30 30 10],'Transform', {'Rotate_Z',alpha});
regards
Thorsten

sebastian
Posts: 114
Joined: Mon 27 Jun 2011, 12:36
Contact:

Re: Oblique emission and rotation troubles

Post by sebastian » Mon 26 Mar 2012, 22:23

Dear Youry,

if I understand correctly, you want to have a position dependent excitation phase. This is currently not possible.
To approximate it, you could define a sinusoidal excitation and create multiple excitation boxes at the correct positions.
For one frequency and medium a phase shift can be transformed into a space shift and vice versa.

A position dependent excitation function should be easy to implement using the extension concept in openEMS (you need to modify the c++ sources). Take engine_ext_excitation.cpp as an example.

kind regards,
Sebastian

Youry
Posts: 51
Joined: Wed 22 Feb 2012, 12:50

Re: Oblic emission and rotation troubles

Post by Youry » Tue 27 Mar 2012, 03:47

Thorsten, first thank you for the missing ' ' in my Rotate_Z, I use to much example without thinking... I should find it be myself, but anyway, it will help me a lot...
Then what you propose to me to do is to put a alpha rotation to every object.
So then each object will rotate around their own gravity center (I imagine), but it is not really what I want. I want that all the object rotate together (like space rotation) without changing anything on the structure.

I feel it can be very very convenient to have 'Transform' possibility on group object for the permittivity geometry definition. At least translation and rotation...
Some boolean operation between object like intersection, union, difference can be also very nice...
But hard work to do I guess...

sebastian, Thank you for your advice, in fact it is exactly what I wanted to do... That why I try to rotate the plan emission... But I think it doesn't work like I want (I will try quickly, now I get correct usage of rotation).
Problem is : If I define a plan (box with 1 dimension 0) and then rotate it, I am scare that very few point (maybe only 1 point, or 1 line) will be define (intersection between grid and arbitrary continuous plan orientation...).

So the technic to define point by point seems good idea, but I am little scare to touch the c++ code. So I will probably do it inside a Octave function...
Last edited by Youry on Tue 27 Mar 2012, 09:08, edited 1 time in total.

Youry
Posts: 51
Joined: Wed 22 Feb 2012, 12:50

Re: Oblic emission and rotation troubles

Post by Youry » Tue 27 Mar 2012, 08:37

So here the point to point phase dependent definition for strange propagation direction.
I define a plane (xy direction) and I put the good phase thank to SetCustomExcite at each point of my plan, to get mainly z propagation, and little bit along x or y...
Still don't work, I don't know why...

I did this :

Most important line are :

Code: Select all

for i=1:max(size(indexX))
for j=1:max(size(indexY))
FDTD = SetCustomExcite(FDTD,f0,[ 'sin(2*pi*' num2str(f0) '*t-' num2str(phase(i,j)) ')' ]);
CSX = AddExcitation(CSX,'excitation',0,pola);
CSX = AddPoint(CSX,'excitation',0,[mesh.x(indexX(i)) mesh.y(indexY(j)) z]); %Ajoute l'exitation dans la grille
end
end
I get an error all the times. This is :
Warning: Unused primitive (type: Point) detected in property: excitation

Total sub-function code :

Code: Select all

function [FDTD,CSX] = SetAddOblicSinus(FDTD,CSX,mesh,unit,index,f0,xm,xM,ym,yM,z,theta,phi,pola00)
%Set a sinus generation, in a plan transverse to Z, with phase space depence in rrespect of theta, phi propagation direction
%First theta rotation of the wavevector with axis x
%Phi is angle between plane zx and vecteur
%The plane generation is between [xm xM] and [ym yM] at the z position of the grid
%Grid has to show a constant xstep and ystep, if not => more energy will be present in the direction of higher density grid location...
%The plane emission is inside a medium of refractive index 'index'
%The polarisation is given in the case of theta=0 and phi=0, and next rotate in respect of angle, [0 1 0] 


%Direction of the waves

M1=[1 0 0;0 cos(theta) sin(theta);0 -sin(theta) cos(theta)];
M2=[cos(phi) 0 sin(phi);0 1 0;-sin(phi) 0 cos(phi)];
M=M2*M1;
u=M*[0 0 1]';
ux=u(1);
uy=u(2);
uz=u(3);
pola=(M*pola00')';
%Find the index x of grid which is inside the plane (in indexX array)
d=0;
for i=1:max(size(mesh.x))
if (mesh.x(i)>=xm)
if (mesh.x(i)<=xM)
d=d+1;
indexX(d)=i;
endif
endif
end
%Find the index y of grid which is inside the plane (in indexY array)
d=0;
for i=1:max(size(mesh.y))
if (mesh.y(i)>=ym)
if (mesh.y(i)<=yM)
d=d+1;
indexY(d)=i;
endif
endif
end
%We check the phase value of all the grid point inside the plan emission
k0=2*pi*f0/299792458*index;
kx=k0*ux;
ky=k0*uy;
for i=1:max(size(indexX))
for j=1:max(size(indexY))
rx=mesh.x(indexX(i))*unit;
ry=mesh.y(indexY(j))*unit;
phase(i,j)=rx*kx+ry*ky;
end
end
%We put sinus and phase at each point of the emission
for i=1:max(size(indexX))
for j=1:max(size(indexY))
FDTD = SetCustomExcite(FDTD,f0,[ 'sin(2*pi*' num2str(f0) '*t-' num2str(phase(i,j)) ')' ]);
CSX = AddExcitation(CSX,'excitation',0,pola);
CSX = AddPoint(CSX,'excitation',0,[mesh.x(indexX(i)) mesh.y(indexY(j)) z]); %Ajoute l'exitation dans la grille
end
end

endfunction
And I call this function with my m-file as follow :

Code: Select all

close all
clear
clc

Xsize=40;
Ysize=40; 
Zsize=40;
unit = 10e-6;

% init and define FDTD parameter
FDTD = InitFDTD(500,1e-5,'CoordSystem',0);  

BC = {'PML_8' 'PML_8' 'PEC' 'PEC' 'PML_8' 'PML_8' }; 
FDTD = SetBoundaryCond(FDTD,BC);

% init and define FDTD mesh
CSX = InitCSX();

mesh.x =0:Xsize;
mesh.y =0:Ysize;
mesh.z =0:Zsize;

CSX = DefineRectGrid(CSX, unit, mesh);  
xm=15;
xM=25;
ym=3;
yM=13;
z=4;
theta=pi/4;
phi=0;
f0=1e12;
pola00=[0 1 0];
[FDTD,CSX]=SetAddOblicSinus(FDTD,CSX,mesh,unit,1,f0,xm,xM,ym,yM,z,theta,phi,pola00);

CSX = AddDump(CSX,'Probe','DumpMode',0);
CSX = AddBox(CSX, 'Probe', 0, [20 0 0],[20 Ysize Zsize]);
rmdir('tmp','s');mkdir('tmp');
WriteOpenEMS('tmp/tmp.xml',FDTD,CSX);
CSXGeomPlot( 'tmp/tmp.xml' );
RunOpenEMS('tmp','tmp.xml','');
disp('use Paraview to visualize the FDTD result...');
Did I do it in bad manner ? Sebastian, can you be little bit more explicit ?

Edit : With AddBox instead of addPoint, I get something, but not what I exepted... Maybe some bug in the function as I wrote quickly and not debug yet... Maybe I will be able to do what I want...

thorsten
Posts: 1438
Joined: Mon 27 Jun 2011, 12:26

Re: Oblic emission and rotation troubles

Post by thorsten » Tue 27 Mar 2012, 11:18

Hi all,

I think I need to clarify a few points.

1) Every transformation is just a plain coordinate transformation. This means, a rotation is always around the point of origin (x=0,y=0,z=0) not on an objects center point or something. Use CSXGeometryPlot to view the effects on your structures.
(If you want it to be around the center of the object, Translate the center to origin first, rotate and translate back)

2) You already can define a position dependent excitation amplitude. Have a look at the Excitation Weighting scheme. e.g. in http://www.openems.de/index.php/Tutoria ... _Waveguide

3) FDTD Excitation Signal means the globally used excitation time-signal:
Most important line are :

Code: Select all

Code: Select all
    for i=1:max(size(indexX))
    for j=1:max(size(indexY))
    FDTD = SetCustomExcite(FDTD,f0,[ 'sin(2*pi*' num2str(f0) '*t-' num2str(phase(i,j)) ')' ]);
    CSX = AddExcitation(CSX,'excitation',0,pola);
    CSX = AddPoint(CSX,'excitation',0,[mesh.x(indexX(i)) mesh.y(indexY(j)) z]); %Ajoute l'exitation dans la grille
    end
    end
This will not work, since there can only be one global time-signal, e.g. a gaussian pulse.
Just use a Sinusoidal Signal as it is predefined...

Code: Select all

CSX = AddExcitation(CSX,'excitation',0,pola,'Delay', 1/f0/4);
This is just where you want to excite (boxes, points e.g.) and with what amplitude and time-delay.

4) Keep in mind the FDTD Yee-Cells scheme... If you excite an electric field, it must be defined on the edges of a cell, not on the nodes ... Therefore a point on the nodes of the mesh will never hit an edge with an electric field to excite... This explains your unused primitives...

I hope I didn't miss anything...

Try using a rotated box (plane) with some thickness (maybe 2-3 cells thick) and excite it uniformly. I don't know if this will work...

It is not implemented yet, but it of course would be useful not only to have a position weighted excitation amplitude, but also time-delay (to mimic phase shift).

Warning: Keep in mind that a time-delay to mimic a phase shift only works for a single frequency... Which is why it is rarely used in FDTD, since the strongest advantage of FDTD is it's broadband excitation capabilities.
If you are interested in a single frequency only, a FD tool (e.g. FEM) would maybe be better suited for your problem...

regards
Thorsten

Youry
Posts: 51
Joined: Wed 22 Feb 2012, 12:50

Re: Oblic emission and rotation troubles

Post by Youry » Tue 27 Mar 2012, 12:20

1) Every transformation is just a plain coordinate transformation. This means, a rotation is always around the point of origin (x=0,y=0,z=0) not on an objects center point or something. Use CSXGeometryPlot to view the effects on your structures.
(If you want it to be around the center of the object, Translate the center to origin first, rotate and translate back)
Ok, I didn't know it, thank you... You right, I should try more, so if I made an alpha transform to every object, all the object will move... That very convenient style to define every rotation around 0 0 0 center... I see.

4) Keep in mind the FDTD Yee-Cells scheme... If you excite an electric field, it must be defined on the edges of a cell, not on the nodes ... Therefore a point on the nodes of the mesh will never hit an edge with an electric field to excite... This explains your unused primitives...
Oh yes that right :oops:
Try using a rotated box (plane) with some thickness (maybe 2-3 cells thick) and excite it uniformly. I don't know if this will work...
Ok I will try this tomorrow... If it failed I will try with 'Delay' option...

This will not work, since there can only be one global time-signal, e.g. a gaussian pulse.
Just use a Sinusoidal Signal as it is predefined...

I don't understand, t is global for every cell, just the phase is different depending on position. If it can work, I will be able to put Fourier series long enough to covert the time calculation duration, and get continuous spectrum...

Anyway, thank you for explanation on rotation, and advice about strange emission, I will have to look this point in more detail soon...

thorsten
Posts: 1438
Joined: Mon 27 Jun 2011, 12:26

Re: Oblic emission and rotation troubles

Post by thorsten » Tue 27 Mar 2012, 16:28

The time signal applied by openEMS (e.g. a Gaussian pulse) is defined just once for the entire simulation. Let's call it e(t).
You can find this signal dumped into the file "et" in the simulation folder and read it just as every other voltage or current definition.

When you define an excitation as a box or some other geometry, you just define a (local) amplitude and optional a delay.
For example you define a vector [1 0 0] (as Volt per meter) as your excitation. This vector is multiplied by e(t - local_Delay) and that's it.

This way, it is ensured that you have the same time-signal (e.g. Gaussian pulse) everywhere. Otherwise it would be a total mess, if you would have a sinusoidal excitation somewhere and a Gaussian pulse somewhere else...

I really need to start working on this User-Manual, I just don't know how to find some time for it... :(

regards
Thorsten

Post Reply