Re: [Meep-discuss] Gaussian Beam source

2020-11-10 Thread Ian Sage


On 09/11/2020 03:54, Ardavan Oskooi wrote:

sources = [mp.GaussianBeamSource(src=mp.ContinuousSource(fcen),
 center=mp.Vector3(),
 size=mp.Vector3(s,s,0),
 beam_x0=beam_x0,
 beam_kdir=beam_kdir,
 beam_w0=beam_w0,
 beam_E0=mp.Vector3(1,1j,0))] 



Many thanks for your help, Ardavan. This works like a charm in 3d. I'd 
like to use 2d for efficiency, and so far my efforts fail. In 2d I have 
a source definition like:


beam_x0 = mp.Vector3(0,0,0)  # beam center
rot_angle = math.radians(0)  # CCW rotation angle about z axis (0: +y axis)
beam_kdir = mp.Vector3(0,1,0).rotate(mp.Vector3(0,0,1),rot_angle) # beam 
propagation direction

beam_w0 = 0.8  # beam waist radius
fcen = 1
df=0.1

sources = [mp.GaussianBeamSource(src=mp.ContinuousSource(frequency=fcen, 
fwidth=df),

 center=mp.Vector3(0,0,0),
 size=mp.Vector3(s,0,0),
 beam_x0=beam_x0,
 beam_kdir=beam_kdir,
 beam_w0=beam_w0,
 beam_E0=mp.Vector3(1,0,0)]

If I set beam_E0=mp.Vector3(1,0,0) or beam_E0=mp.Vector3(0,0,1) it works 
fine, but if E0 is set to any other direction such as 
beam_E0=mp.Vector3(1,0,1) or beam_E0=mp.Vector3(1j,0,1) then I end up 
with fields of zero. Normalizing |E0| to 1 doesn't help.


Is there any way to generate a circular polarized Gaussian beam in 2d? 
(it's possible with a simulated QWP, of course)


Many thanks,

Ian


___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

[Meep-discuss] Tiny bugs in build script (I think)

2020-11-09 Thread Ian Sage
I just needed to do a fresh Meep install, on a fresh install of Ubuntu 
18.04. As usual, I used the script in the "build from source" section of 
the docs. All was well, except:


1) The meep python package is installed to 
/usr/local/lib/python3.6/site-packages but that directory is not added 
to the PYTHONPATH environment variable. I guess for many users it will 
already have been added by installation of other packages.


2) I had to retain the line:
    sudo pip3 install --upgrade pip
which is marked as only needed for Ubuntu 16.04. If that line is 
commented out the build completes successfully but it leads to "meep has 
no Vector3 object" errors.


Hope that helps someone,

Ian


___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

[Meep-discuss] Gaussian Beam source

2020-11-08 Thread Ian Sage
I was prompted by Steven Johnson's recent posting (1st Nov) to switch
some of my simulations to use the MEEP GaussianBeamSource. The
documentation states that a complex valued polarization vector can be
used to generate a circular polarized beam, but there are no explicit
instructions for this. I seem to have tried every sane combination of y
and z values (for an x-axis propagating beam, source is a line along the
y-axis, 2-d simulation) without success.

Does anyone have a working sample circular polarized GaussianBeamSource
code, which they can share?

Many thanks,

Ian


___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

Re: [Meep-discuss] GaussianBeamSource

2020-10-31 Thread Ian Sage
You just need to apply an amp_func function to your source, with an
appropriate spatial profile. The function takes as argument a vector
giving the spatial offset from the centre of the source, allowing you to
specify a scaled exponential fall-off from the centre. You can also
return a complex amplitude to spatially vary the phase and steer the beam.

It's pretty well covered in the documentation.

Ian

On 31/10/2020 16:20, Alfredo Daniel Sánchez wrote:
> Hello everybody!
>
> I realized that the GaussianBeamSource class is not defined in my meep
> package. It is possible that it does'n exist anymore? In that case,
> how can I generate one?
>
> Best regards,
>
> Alfredo.
>
> ___
> meep-discuss mailing list
> meep-discuss@ab-initio.mit.edu
> http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

Re: [Meep-discuss] Providing geometry as voxel-array via python

2019-02-21 Thread Ian Sage
You can use a material function to do this from Python. - 
straightforwardly if your voxel data matches the resolution of the Meep 
simulation but you'd need to interpolate or discretize otherwise


Ian

On 21/02/2019 14:13, Jonathan wrote:

Hello,

the normal way to set up the scene geometry via the Python interface 
is to build it out of primitives (boxes, spheres, etc.). I however 
already have the permitivities as a 3d array. The python interface 
does not seem to provide a direct way to set them. I think it is maybe 
possible via the low level C++ interface, however I could not find any 
examples for this. What is the recommended way (considering 
robustness, performance, etc...) to do so?


With kind regards,

Jonathan Klein

___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss


___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

[Meep-discuss] Class methods as materials functions

2019-01-29 Thread Ian Sage
I am able to tailor Meep sources by using a src_func, and have no 
problem writing that function either as a plain standalone, or as a 
class method. The latter is often handy if only to avoid a proliferation 
of global variables etc. Equally, I can create custom materials by 
providing a material_function - provided I write this as a standalone. 
If I try to implement a material_function as a class method, I see an 
error thrown by the line defining the GeometricObject, which ends in:


File "/usr/local/lib/python3.5/site-packages/meep/geom.py", line 391, in 
__init__

    super(Block, self).__init__(**kwargs)
  File "/usr/local/lib/python3.5/site-packages/meep/geom.py", line 306, 
in __init__

    material.eps = False
AttributeError: 'method' object has no attribute 'eps'

Is it possible to use a class method as a materials_function? Do I need 
a different idiom?


Many thanks for any advice,

Ian


___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

Re: [Meep-discuss] Problem using epsilon_offdiag with load_minus_flux_data

2018-12-21 Thread Ian Sage
On 21/12/2018 18:56, Ardavan Oskooi wrote:
> You will need to add "force_all_components=True"

Thank you, Ardavan - that has fixed it.

As a matter of interest, can I now use that flag in a 1-d simulation to
access E_y field components?

Many thanks,

Ian


___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

[Meep-discuss] Problem using epsilon_offdiag with load_minus_flux_data

2018-12-21 Thread Ian Sage
I have been trying to examine the reflection of light from a
(semi-infinite) birefringent slab, and encounter a problem when the
optic axis is rotated, giving non-zero entries in epsilon_offdiag. Under
this condition, the program fails at the load_minus_flux_data command,
with the error message:

    meep: Total dft_chunks size does not agree with size allocated for
output array.

emitted by each thread. The fairly minimal example attached, shows this
problem.

The error can be cleared by un-commenting line 29 (to set the optic axis
rotation to zero) or line 39 (to artificially set epsilon_offdiag to
zeros) or by commenting out line 70 (to avoid calling
load_minus_flux_data). In the latter case, sensible results are given by
performing the flux subtraction manually, outside meep.

Am I doing something wrong, or is this a tiny bug?

Many thanks for any help, and best wishes to all for the holiday,

Ian







# Model reflection at normal incidence
import math
import meep as mp
import numpy


layer_thickness=25.0
ne=1.7
no=1.5

lambda_min=0.4
lambda_max=0.6
f_min=1/lambda_max
f_max=1/lambda_min
freq=(f_min+f_max)/2
df=f_max-f_min

sim_x=1
sim_y=layer_thickness
dpml=1
res=25

fn_call_count=0

def mat_function(pos):
global fn_call_count
fn_call_count=fn_call_count+1
rot_angle=math.pi/4
#rot_angle=0.0
c=math.cos(rot_angle)
s=math.sin(rot_angle)
Q=numpy.matrix([[c, 0, -s],[0, 1, 0],[s, 0, c]]) # Rotation matrix around y
Q_inv=Q.transpose()
T_eps=numpy.matrix([[ne*ne, 0, 0],[0, no*no, 0],[0, 0, no*no]])
p1=T_eps*Q_inv
p2=Q*p1
diag=mp.Vector3(p2[0, 0], p2[1, 1], p2[2, 2])
off_diag=mp.Vector3(p2[0, 1], p2[0, 2], p2[1, 2])
#off_diag=mp.Vector3(0.0, 0.0, 0.0)
if fn_call_count%1==0: # Print out just a few values from the run
print(pos, '\n', rot_angle, '\n', Q, '\n\n', Q_inv, '\n\n', T_eps, '\n\n', p1, '\n\n', p2, '\n\n', diag, '\n', off_diag, '\n**\n')
return mp.Medium(epsilon_diag=diag, epsilon_offdiag=off_diag)


cell=mp.Vector3(sim_x, sim_y+2*dpml, 0)
pml=mp.PML(dpml, direction=mp.Y)


src=mp.Source(src=mp.GaussianSource(frequency=freq, fwidth=df), component=mp.Ex, size=mp.Vector3(cell.x, 0, 0), center=mp.Vector3(0, -0.5*sim_y, 0))
geom=[]

file_suffix='reference'
sim=mp.Simulation(cell, res, geometry=geom, sources=[src], boundary_layers=[pml], dimensions=2, k_point=mp.Vector3(0, 2, 0))

detector=mp.FluxRegion(center=mp.Vector3(0, -0.5*(sim_y-1), 0), size=mp.Vector3(sim_x, 0, 0))
ref_flux=sim.add_flux(freq, df, 101, detector)

sim.run(mp.to_appended(file_suffix, mp.at_every(0.5, mp.output_efield_x)), until_after_sources=2*sim_y)
source_ref=sim.get_flux_data(ref_flux)
sim.display_fluxes(ref_flux)

sim.reset_meep()
geom=[]
biref_layer=mp.Block(size=mp.Vector3(cell.x, layer_thickness, 1e20), center=mp.Vector3(0, 1, 0), material=mat_function)
geom.append(biref_layer)
file_suffix='main_run'
sim=mp.Simulation(cell, res, geometry=geom, sources=[src], boundary_layers=[pml], dimensions=2, k_point=mp.Vector3(0, 2, 0))
detector2=mp.FluxRegion(center=mp.Vector3(0, -0.5*(sim_y-1), 0), size=mp.Vector3(sim_x, 0, 0))
whole_flux=sim.add_flux(freq, df, 101, detector2)
sim.load_minus_flux_data(whole_flux, source_ref)
sim.run(mp.to_appended(file_suffix, mp.at_every(0.5, mp.output_efield_x)), until_after_sources=2*sim_y)
sim.display_fluxes(whole_flux)
___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

[Meep-discuss] Chiral LC simulation - a question of dimension

2018-08-29 Thread Ian Sage
As one of my main interests is in liquid crystal layer optics, I have 
recently tried to model the classic textbook example of reflection from 
a cholesteric material. Such a material is optically uniaxial, with a 
helical arrangement of the optic axis, with a pitch length in the 
optical range. I expect to see 50% reflection of linear polarized 
incident light, within a narrow band centred on the frequency 
corresponding to wavelength in the LC equal to the helical pitch length. 
(The reflection is circular polarized.)


If I carry out a 2D or 3D simulation, this is indeed what I observe - 
although there is a strong and strange wavelength dependency 
superimposed, which I need to check through. But if I use a 1D 
simulation, I see a 100% reflection coefficient (and also lose the 
anomalous wavelength effect).


Is there some issue with 1D simulations which means they cannot cope 
with a helical permittivity tensor?


Many thanks for any guidance available,

Ian

___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

[Meep-discuss] Install from git clone - another fail

2018-08-07 Thread Ian Sage
On Arch linux, installing from git through the AUR packages, I now see a
problem reminiscent of (but different from) John's:

In file included from meep-python.cpp:3794:
typemap_utils.cpp: In function ‘std::__cxx11::string
py_class_name_as_string(PyObject*)’:
typemap_utils.cpp:21:51: error: invalid conversion from ‘const char*’ to
‘char*’ [-fpermissive]
 #define PyObject_ToCharPtr(n) PyUnicode_AsUTF8(n)
   ^~~
This follows an update of harminv, libctl, h5utils, mpb without
problems, motivated by a compiler libraries update.

Any clues?

Many thanks,

Ian

On 02/08/18 06:24, John Weiner wrote:
> Meep.hpp:1275:78 error: invalid conversion from ‘meep::component*’ to
> ‘int’ -fpermissive:
> chunk loop filed components (fc, grid, shift, phase, S, sh,
> components.data(), components.size()) {}
>
>

___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

Re: [Meep-discuss] Neophyte looking for example of continuous symmetry.

2018-05-08 Thread Ian Sage
I am learning too! There is a write-up of the difference between
periodic and Bloch periodic boundaries at

https://kb.lumerical.com/en/ref_sim_obj_bloch_bc.html

which I find more helpful than other accounts.

Ian


On 09/05/18 00:08, Doug McKnight wrote:
> Hi Joe,
> I think I'll need to spend some time experimenting at this point.
> Thanks,
> Doug
>
> On 5/8/2018 2:35 PM, Joe Lowney wrote:
>> Hi Doug,
>>
>> Your structure will be, by definition, periodic. The periodicity will
>> be whatever the spatial extent of the geometry that you defined is.
>>
>> Your fields will also be periodic, by definition, unless you specify
>> bloch periodic boundary conditions - which introduce a phase offset
>> at the boundary. If you set the k-point to be the k-vector of your
>> plane wave, this phase offset will match the phase in the way you want.
>>
>> In other words, you don't have to do anything to make your structure
>> repeating: this is automatic.  You also don't have to tell MEEP that
>> this periodicity goes on forever, this is already true.   You
>> /do/ need to implement bloch-periodic boundary conditions in order to
>> phase match the fields for an arbitrary input wave.  If you set the
>> k-point to be the input plane wave k-vector, you should be good to go
>> /for any input angle/.
>>
>> You do not need to manually choose the angles to get the phase to
>> match: this is the advantage of Bloch periodic boundary conditions.
>>
>> Hope this helps.
>>
>> Joe
>>
>>
>>
>>
>>
>> On Tue, May 8, 2018 at 12:43 PM Doug McKnight <d...@mcknight.to
>> <mailto:d...@mcknight.to>> wrote:
>>
>> Joe,
>> Thanks for the reply.  Maybe this is the source of my confusion:
>> In the docs, the k_point vector defines the periodicity of the
>> *fields* at the boundaries.  But, I want to have a repeating
>> structure with a, possibly different, periodicity.
>>
>> For instance (2-D case), if I make a thin diffraction grating
>> which is repeats every /n/ units in y then, intuitively, it seems
>> I should make my computation cell of size /n/ units in the y
>> direction, and somehow tell Meep that it repeats forever...
>>
>> Then, if I illuminate it with a plane wave at some arbitrary
>> angle, the Bloch condition isn't necessarily satisfied, right? 
>>
>> Is this what you mean by forcing "true periodic boundary
>> conditions", I need to manually choose the angles to get the
>> phase to match?
>>
>> Apologies if I'm being dense.
>>
>>
>> Doug
>>
>>
>>
>>
>> On 5/8/2018 12:54 PM, Joe Lowney wrote:
>>> Doug - 
>>>
>>> If you've defined a general plane wave (kx ky kz), then the
>>> bloch vector, k-point, can be set to (kx ky kz) / (2 * pi). 
>>> This should automatically handle phase matching at the
>>> boundaries for an arbitrary angle of incidence.  You shouldn't
>>> have to replicate any unit cells or anything complicated.
>>>
>>> The discrete angle choices would be required if you were trying
>>> to force-match true periodic boundary conditions, but this
>>> shouldn't be necessary with bloch-periodic BCs, which allows
>>> your incident angle to be completely arbitrary.
>>>
>>> Joe
>>>
>>>
>>>
>>>
>>> On Tue, May 8, 2018 at 11:41 AM Doug McKnight <d...@mcknight.to
>>> <mailto:d...@mcknight.to>> wrote:
>>>
>>> Ian and Joe,
>>> Thanks, yes, that's the sort of approach I'm attempting to
>>> take. It's
>>> just a little slow while I find my way into the Meep way of
>>> thinking... 
>>> It would probably help if I was a "real programmer".
>>>
>>> I was presuming that, at least initially, I'd need to
>>> phase-match my
>>> incident wave at the boundary. To achieve a more fine-tuned
>>> choice of
>>> incident angles I could replicate n copies of my unit cell
>>> within the
>>> bloch-periodic simulation region. Then the phase-matching
>>> condition
>>> would still force discrete angle choices, but they'd be more
>>> closely
>>> spaced.
>>>
>>> If there's a more elegant approach, I'm all ears.
>>> Tha

Re: [Meep-discuss] Neophyte looking for example of continuous symmetry.

2018-05-08 Thread Ian Sage
> so far, I'm happily sending an oblique plane wave at a block of 
glass, and seeing it refract/reflect as expected.


> I'd like to transform my finite block of stuff to be infinitely long 
and specified in terms of a unit cell, (and the same for the
> plane-wave source) so that I can construct diffraction-grating-like 
features.


So, you want to use an oblique source with cyclic boundary conditions? 
Are you working in 2d or 3d?


Either way, you'll need to remove the PML from the axis/axes along which 
you want periodicity, using the direction parameter, and specify a 
k_point vector. Both points are covered in the Python UI docs.


If you want off-axis incidence, you'll also need to adjust the 
wavelength/simulation dimension/incidence angle to ensure smooth phase 
matching across the cyclic boundary (unless a real expert knows a better 
way).


Is that what you wanted to know?

Ian


___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

[Meep-discuss] Python material_function and extra_materials

2018-04-09 Thread Ian Sage
I have no problems using a material_function to set an index which
varies with position. However, the construct fails if my function
returns materials_library.Al even if I include this in an
extra_materials list (either in the initialization of the simulation or
as a separate assignment). The simulation runs without error, but the
metal is seemingly absent.

All is well if I just assign a Block material as materials_library.Al

Are these (material_function; extra_materials) constructs fully
implemented? If so, do they need a special syntax?

I am running meep-git-1.4.3.r76 on Arch.

Thanks,

Ian


___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

Re: [Meep-discuss] Custom Amp Function for Chirping

2018-03-28 Thread Ian Sage
On 28/03/18 02:21, Priscilla Kelly wrote:
> Hello,
>
> I am trying to add a phase in terms of time but by setting the amplitude as a 
> function. 
>
> I see there is a way to do it as a function of space, but how do you do it as 
> a function of time? I’d like to add a chirp to the pulse.
I cannot help you with Scheme implementation, but it is easy in Python.

Ian


import meep as mp
import math
import cmath

sim_x=5
sim_y=5
dpml=1.0

res=50

t=150

start_f=0.5
end_f=2
duration=100

def source_function(t):
f=start_f+(end_f-start_f)*t/duration
return cmath.exp(complex(1, -2*math.pi*f*t))

src=[mp.Source(mp.CustomSource(src_func=source_function, end_time=duration), component=mp.Ez, center=mp.Vector3(0, -0.5*sim_y, 0), size=mp.Vector3(sim_x+2*dpml, 0, 0))]
pmls=[mp.PML(dpml)]
cell=mp.Vector3(sim_x+2*dpml, sim_y+2*dpml, 0)

sim=mp.Simulation(cell_size=cell, sources=src, boundary_layers=pmls, resolution=res, dimensions=2, force_complex_fields=True)
sim.run(mp.to_appended('ez', mp.at_every(0.1, mp.output_efield_z)), until=t)
___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

[Meep-discuss] Documentation of use_output_directory

2018-03-25 Thread Ian Sage
In the Python Tutorial/Basics documentation, it says:

Instead, we can add the following command to run:
mp.use_output_directory()

but my attempts to use the command in this way (as a clause in the
sim.run command) have been unsuccessful. Everything is fine if I do:

sim.use_output_directory()

outside the run command. The account of the command in the Python user
interface section is closer to my experience, but I still read it as
suggesting that the command is an attribute of meep, rather than
meep.Simulation.

HTH,

Ian


___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

Re: [Meep-discuss] Argument required for output_efield()?

2018-03-24 Thread Ian Sage
Ah, OK -

sim.run(mp.with_prefix('fields', mp.in_point(mp.Vector3(0, 0, 0),
mp.at_time(24, mp.synchronized_magnetic(mp.output_efield,
mp.output_hfield, until=t)

seems to work fine. There is an apparent anomaly in that using
meep.with_prefix as above, puts the e-field and h-field data in separate
files, while if I substitute a meep.to_appended in its place, all the
components are output to the same file.

As an aside, I am probably not the first to be confused by the
presentation in the documentation for the Python user interface*||*|,|
of the output functions (*|output_epsilon() etc)|*, |as typical function
calls||when ||- as is clear from the tutorials - a function
reference is what is needed.|*|
|*

*|Ian
|*


On 23/03/18 23:54, Ian Sage wrote:
> I am trying to analyze the polarization state of light at a particular
> point in a simulation. As the propagation direction may vary, I plan to
> record the three components of both the electric and magnetic fields as
> a first step. I guess I need to apply synchronize_magnetic, at_time and
> in_point to specify the point and get accuracy? So I have been trying to
> build up the run/step function:
>
> sim.run(mp.synchronized_magnetic(mp.in_point(mp.Vector3(0, 0, 0),
> mp.at_time(24, mp.with_prefix('le', mp.output_efield().
>
> however, I get an error thrown:
>
> TypeError: output_efield() missing 1 required positional argument: 'sim'
>
> The documentation does not seem to specify any argument to
> output_efield() - what construct should I be using?
>
> More generally, is there a better solution to outputting all the field
> components at a specific time and point?
>
> Many thanks,
>
> Ian
>
>

___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

[Meep-discuss] Argument required for output_efield()?

2018-03-23 Thread Ian Sage
I am trying to analyze the polarization state of light at a particular
point in a simulation. As the propagation direction may vary, I plan to
record the three components of both the electric and magnetic fields as
a first step. I guess I need to apply synchronize_magnetic, at_time and
in_point to specify the point and get accuracy? So I have been trying to
build up the run/step function:

sim.run(mp.synchronized_magnetic(mp.in_point(mp.Vector3(0, 0, 0),
mp.at_time(24, mp.with_prefix('le', mp.output_efield().

however, I get an error thrown:

TypeError: output_efield() missing 1 required positional argument: 'sim'

The documentation does not seem to specify any argument to
output_efield() - what construct should I be using?

More generally, is there a better solution to outputting all the field
components at a specific time and point?

Many thanks,

Ian


___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

[Meep-discuss] Possible bug in output_png parallel operation?

2018-03-10 Thread Ian Sage
I am running a simple 2d simulation using the latest de8571e commit code
on Arch Linux. The python file is named gaussian_beam.py and output is
by the output_png instruction:

sim.run(mp.at_beginning(mp.output_epsilon), mp.at_every(0.1,
mp.output_png(mp.Ez, '-c bluered -Z')), until=t)

If I do: "python gaussian_beam.py" or "mpirun -N 1 python
gaussian_beam.py", all is well. But if I use "mpirun -N 2 python
gaussian_beam.py", I see a huge sequence of errors, beginning:

creating output file "./gaussian_beam-eps-00.00.h5"...
creating output file "./gaussian_beam-ez-00.10.h5"...
HDF5-DIAG: Error detected in HDF5 (1.10.1) thread 0:
  #000: H5F.c line 586 in H5Fopen(): unable to open file
major: File accessibilty
minor: Unable to open file
  #001: H5Fint.c line 1305 in H5F_open(): unable to lock the file
major: File accessibilty
minor: Unable to open file
  #002: H5FD.c line 1839 in H5FD_lock(): driver lock request failed
major: Virtual File Layer
minor: Can't update object
  #003: H5FDsec2.c line 940 in H5FD_sec2_lock(): unable to lock file,
errno = 11, error message = 'Resource temporarily unavailable'
major: File accessibilty
minor: Bad file ID accessed
h5topng error: error opening HD5 file
creating output file "./gaussian_beam-ez-00.20.h5"...
HDF5-DIAG: Error detected in HDF5 (1.10.1) thread 0:

All the error messages seem to come from thread 0.
If I change the step function to write to hdf5:

sim.run(mp.at_beginning(mp.output_epsilon), mp.to_appended("ez",
mp.at_every(0.1, mp.output_efield_z)), until=t)

there is no error reported.

I am not certain, whether this is a new issue with the latest commit, or
if it also affects earlier versions.

Hope that's useful to someone!

Ian

___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

[Meep-discuss] Problem using output_field_function in at_every

2018-02-07 Thread Ian Sage
The documentation is not all that clear about the procedure to follow,
when using the Python interface. I'd like to call a function like:

def esq_mag(r, ex, ey, ez):
    return ex*ex.conjugate()+ey*ey.conjugate()+ez*ez.conjugate()

from the at_every construct - I initially used the above function
together with:

sim.run(mp.to_appended(prefix, mp.at_every(0.1,
sim.output_field_function('E_sq', [mp.Ex, mp.Ey, mp.Ez], esq_mag))),
until=t)

but the program runs to the point where output is triggered, then stops
with a lengthy error message, terminating:

      File "/usr/lib/python3.6/site-packages/meep/simulation.py", line
28, in get_num_args
        return func.__code__.co_argcount
    AttributeError: 'NoneType' object has no attribute '__code__'

That does not surprise me as the docs state I have to pass a function to
the sim.run function, not the result of evaluating a function. So I
tried again, supplying in the spirit of the instructions related to
Scheme "a function |my-weird-output| of no arguments that, when called,
outputs our function |f|." So I constructed a function to do that:

def esq_mag():
    def _esq_mag(r, ex, ey, ez):
    return ex*ex.conjugate()+ey*ey.conjugate()+ez*ez.conjugate()
    return _esq_mag

but the result is another error:

    TypeError: esq_mag() takes 0 positional arguments but 4 were given

OK - I need to handle those arguments to esq_mag, so

def esq_mag(r, ex, ey, ez):
    def _esq_mag(r, ex, ey, ez):
    return ex*ex.conjugate()+ey*ey.conjugate()+ez*ez.conjugate()
    return _esq_mag

but now the program fails with:

    TypeError: must be real number, not function

    The above exception was the direct cause of the following exception:

    SystemError: PyEval_EvalFrameEx returned a result with an error set

I have tried several variations without success; please, what construct
should I be using?

Also, if I try to use output_real_field_function as described in the
docs, I get the error:

    AttributeError: 'Simulation' object has no attribute
'output_real_field_function'

I'll be very grateful for any help,

Ian

___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

[Meep-discuss] Build fail using latest commit

2018-02-02 Thread Ian Sage
FYI, in my hands using Arch Linux, build of the latest (1st Feb) commit
failed with the following message:

Making all in python
make[2]: Entering directory '/localpath/meep-git/src/meep/python'
swig -Wextra -I../src -I../libmeepgeom -I..  -outdir . -c++ -python -o
meep-python.cpp ./meep.i
cp meep.py __init__.py
if [[ "3.0.12" = 3.0.12 ]]; then \
    sed -i "" '/^if _swig_python_version_info >= (2, 7, 0):/,/^else:/d'
__init__.py; \
    sed -i "" 's/    import _meep/from . import _meep/' __init__.py; \
fi
sed: can't read /^if _swig_python_version_info >= (2, 7, 0):/,/^else:/d:
No such file or directory
sed: can't read s/    import _meep/from . import _meep/: No such file or
directory
make[2]: *** [Makefile:1154: __init__.py] Error 2
make[2]: Leaving directory '/localpath/meep-git/src/meep/python'
make[1]: *** [Makefile:505: all-recursive] Error 1
make[1]: Leaving directory '/local/meep-git/src/meep'
make: *** [Makefile:414: all] Error 2
==> ERROR: A failure occurred in build().
    Aborting...

Hope that is useful.

Ian


___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

Re: [Meep-discuss] Defining Gaussian beam in pymeep

2018-01-20 Thread Ian Sage
> I have tried the following for a beam propagating along the x-axis (so
> y-coordinate should follow a Gaussian profile):

> But it seems to give something more like a point source rather than a
beam.

You don't give all the parameters you used; if the source area is
point-like, it will diverge like any point source (in the attached, set
radius to eg 0.1). You need a substantial area.

Depending on your application you might also want to control the phase
to put some curvature in the wavefront at the source.

Ian
import meep as mp
import math
import cmath

x_size=30
y_size=30

cell=mp.Vector3(x_size, y_size, 0)

dpml=1.0

pml_layers=[mp.PML(dpml)]

resolution=20

radius=2
amp_factor=1
phi_factor=1/4
skew=5.0

def term1(position):
x=position.x
val=amp_factor*(x/radius)**2
return math.exp(-val)

def term2(position):
x=position.x
val=complex(0, x*skew+phi_factor*(x/radius)**2)
return cmath.exp(-val)

def pw_amp():
def _pw_amp(position):
return term1(position)*term2(position)
return _pw_amp

fcen=2.0
df=0.02

sources = [
mp.Source(
mp.ContinuousSource(fcen, fwidth=df),
component=mp.Ez,
center=mp.Vector3(0, -14),
size=mp.Vector3(x_size, 0, 0),
amp_func=pw_amp()
)
]

sim = mp.Simulation(
cell_size=cell,
sources=sources,
boundary_layers=pml_layers,
resolution=resolution,
)

t=50
sim.run(mp.to_appended("ez", mp.at_every(0.1, mp.output_efield_z)), until=t)
___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

[Meep-discuss] PyMeep - output_Xfield_r, output_Xfield_p not implemented?

2018-01-15 Thread Ian Sage
Again in my hands, I get no output or error message when I try to invoke
the above output functions. The corresponding functions
output_Xfield_x,y,z work fine, as expected.

Are the polar functions awaiting implementation?

Thanks,

Ian


___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

Re: [Meep-discuss] PyMeep - possible progress_interval bug?

2018-01-15 Thread Ian Sage
Thanks, Chris.

Ian


On 15/01/18 15:00, Chris Hogan wrote:
> Hi Ian,
>
> Thanks for catching this. I've fixed it in PR#170. 
>
>


___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

Re: [Meep-discuss] Plot field on a line at every time step

2018-01-15 Thread Ian Sage
On 15/01/18 10:22, Kapidani Bernard wrote:

> What is the correct procedure to output the the fields on a physical
> line with a given resolution at each time step?
I cannot say it's a "right" way - and I use the Python rather than
Scheme interface, but in a homebrew step function, looping over the
points with a call to sim.get_field_point seems fine.

Ian

___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

[Meep-discuss] PyMeep - possible progress_interval bug?

2018-01-14 Thread Ian Sage
In my hands, using the 5th Jan development code commit on Arch Linux, an
error:

TypeError: __init__() got an unexpected keyword argument 'progress_interval'

results, if I include an assignment to progress_interval in the
simulation definition:

sim=mp.Simulation(cell_size=cell,
  sources=[source1, source2],
  boundary_layers=pml,
  resolution=res,
  force_complex_fields=True,
  progress_interval=100
)

All is well if I comment out that line & the preceding comma.

Is this a bug or does it need a different syntax?

Thanks,

Ian


___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

[Meep-discuss] Thin metals - further observations

2018-01-10 Thread Ian Sage
I have been looking a little further at the issues surrounding
limitingly thin metal layers. Setting a metal layer thickness to 0.0,
(air/vacuum on each side) still gives a finite reflection. This holds
whether the metal is patterned or continuous across the simulation
volume. It also - apparently - holds very generally whenever any Drude
or Lorentz type term is included in the material definition.

I'd be grateful for any input here - is this just the wrong approach to
use for (finite) thin metals, or a bug? Is there a workaround?

Thanks,

Ian


___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

[Meep-discuss] Simulation corner case - zero thickness metal layer

2018-01-03 Thread Ian Sage
I realize this is abuse of software, but here goes.

I've been attempting to simulate WGP structures and looked at the
transmission as a function of (rectangular) wire thickness at constant
pitch of 200nm (100nm each of wire and gap). If the metal thickness is
taken to zero, the transmission does not go to 1.0 as I'd expect and the
field can be seen interacting with the zero-thickness metal in the
attachment. Seeing these effects for a metal thickness set to zero, I
lose confidence in results using thin, finite metal layers (say, 20nm)

Changing the resolution or turning eps_averaging on/off seems not to
change the outcome. Use of a non-dispersive dielectric material in place
of the Drude Lorentz metal gives the expected, simple behaviour as the
thickness approaches zero.

Is there a "proper" way to include thin metal structures?

Thanks for any advice,

Ian

import sys
import os
#import shutil
from math import pi
import datetime
import meep as mp

susc_aluminium = [
mp.LorentzianSusceptibility(frequency=1e-20, gamma=0.037908, sigma=7.6347e41),
mp.LorentzianSusceptibility(frequency=0.13066, gamma=0.26858, sigma=1941),
mp.LorentzianSusceptibility(frequency=1.2453, gamma=0.25165, sigma=4.7065),
mp.LorentzianSusceptibility(frequency=1.4583, gamma=1.0897, sigma=11.396),
mp.LorentzianSusceptibility(frequency=2.8012, gamma=2.7278, sigma=0.55813)
]

aluminium = mp.Medium(epsilon=1.0, E_susceptibilities=susc_aluminium)

dummy = mp.Medium(epsilon=100.0)


x_size=1.0
y_size=4.0
cell=mp.Vector3(x_size, y_size, 0)
pml_layers=[mp.PML(1.0, direction=mp.Y)]
f_cen=1.9568
df=1.3495
par_source=[mp.Source(mp.GaussianSource(f_cen, fwidth=df), component=mp.Ez, center=mp.Vector3(0, -0.75), size=mp.Vector3(x_size, 0))]
perp_source=[mp.Source(mp.GaussianSource(f_cen, fwidth=df), component=mp.Hz, center=mp.Vector3(0, -0.75), size=mp.Vector3(x_size, 0))]
grating_pitch=0.2
wire_width=0.1
wire_thickness=0.0
grating_geometry=[]
for wire_count in range(5):
grating_geometry.append(mp.Block(mp.Vector3(wire_width, wire_thickness, 1e20), center=mp.Vector3(-0.5*x_size+(0.5+wire_count)*grating_pitch, 0, 0), material=aluminium))
trans_fr = mp.FluxRegion(center=mp.Vector3(0, 0.75), size=mp.Vector3(x_size, 0))
output_dir='WGP-'+datetime.datetime.now().strftime("%Y-%m-%d-%H-%M-%S")
#os.mkdir(output_dir)
trans_fr = mp.FluxRegion(center=mp.Vector3(0, 0.75), size=mp.Vector3(x_size, 0))
sim=mp.Simulation(cell_size=cell, boundary_layers=pml_layers, k_point=mp.Vector3(2*grating_pitch/(2*pi), 0, 0), sources=perp_source, geometry=[], eps_averaging=True, resolution=100)
sim.use_output_directory(output_dir)
tr_flux = sim.add_flux(f_cen, df, 201, trans_fr)
sim.run(mp.at_beginning(mp.output_epsilon), mp.to_appended("hz0", mp.at_every(0.1, mp.output_hfield_z)), until=200)
stdout=sys.stdout
opf=open(output_dir+'/reference_flux.out', 'w')
sys.stdout=opf
sim.display_fluxes(tr_flux)
opf.close()
sys.stdout=stdout
sim.reset_meep()
sim=mp.Simulation(cell_size=cell, boundary_layers=pml_layers, k_point=mp.Vector3(2*grating_pitch/(2*pi), 0, 0), sources=perp_source, geometry=grating_geometry, eps_averaging=True, resolution=100)
sim.use_output_directory(output_dir)
tr_flux = sim.add_flux(f_cen, df, 201, trans_fr)
sim.run(mp.at_beginning(mp.output_epsilon), mp.to_appended("hz", mp.at_every(0.1, mp.output_hfield_z)), until=200)
stdout=sys.stdout
opf=open(output_dir+'/flux.out', 'w')
sys.stdout=opf
sim.display_fluxes(tr_flux)
opf.close()
sys.stdout=stdout

___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

Re: [Meep-discuss] Writing fluxes to file in PyMeep

2018-01-01 Thread Ian Sage

On 31/12/2017 23:45, Ardavan Oskooi wrote:
The save_flux routine outputs Meep's internal data which unfortunately 
is not easily decipherable. We will soon add support for outputting 
the Fourier-transformed fields to a file in a readable format.


Thanks, Ardavan - I'll be on the edge of my seat, waiting for this!
Best wishes for the New Year,

Ian

___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

[Meep-discuss] Writing fluxes to file in PyMeep

2017-12-29 Thread Ian Sage
The display_fluxes() function is very convenient for listing
transmission/reflection spectra, as is well illustrated in the tutorials.

Is there a similar, easy route to write the same information to a file?
The save_flux() function - if I understand correctly - writes out the
Fourier transformed flux information which is not what I need.

I can achieve my aim by redirecting stdout, but is there an easy, direct
way?

Many thanks for any advice,

Ian


___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

Re: [Meep-discuss] Parallel PyMeep - odd effect of resolution change

2017-12-19 Thread Ian Sage

Thanks for your efforts on all our behalves, Ardavan.

Ian


On 19/12/2017 16:34, Ardavan Oskooi wrote:


This is a bug; the patch is on GitHub 
.




___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

[Meep-discuss] Parallel PyMeep - odd effect of resolution change

2017-12-18 Thread Ian Sage
In trying to educate myself about Meep and FDTD, I have been running
various toy simulations, and most recently have tried to replicate metal
reflectivity results as presented at

http://juluribk.com/2011/04/27/plasmonic-materials-in-meep/ ;

the project file linked from that page is also the source of the LD
model parameters I have used. I am using the github source of Dec 8th
commits, on Arch Linux.

According to the value of resolution I choose, the simulation may
succeed or fail with a message such as: "meep: Lorentzian pole at too
high a frequency 1e-20 for stability with dt = 0.00125313: reduce the
Courant factor, increase the resolution, or use a different dielectric
model" but the pattern of success/failure bewilders me. For example:

resolution (set in run_parameters.py)
50 - succeeds
100 - fails
150 - succeeds
200 - succeeds
300 - fails
350 - succeeds
400 - fails
401 - succeeds
399 - fails
398 - succeeds

When the simulation succeeds, the results closely match those in the
source reference and are almost superimposable apart from the
resolution=50 run. I think the simulation files are pretty simple -
attached.

Am I doing something fundamentally wrong - or is there a rational
explanation, or bug?

Many thanks,

Ian



metal_mirror.sh
Description: application/shellscript
import meep as mp
import metals
from run_parameters import *

#k_point=mp.Vector3(0,0,0)
#sym = mp.Mirror(direction=mp.Y, phase=1)

cell=mp.Vector3(x_size, y_size)

sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df),
 component=mp.Ez,
 center=mp.Vector3(0, source_y),
 size=mp.Vector3(x_size, 0))]

pml_layers=[mp.PML(pmld)]

sim=mp.Simulation(cell_size=cell,
  boundary_layers=pml_layers,
  #geometry=geometry,
  sources=sources,
  #symmetries=[sym],
  resolution=resolution)

refl_fr = mp.FluxRegion(center=mp.Vector3(0, fr_y), size=mp.Vector3(x_size/2-pmld, 0))

#pt=mp.Vector3(0, 0)

refl = sim.add_flux(fcen, df, nfreq, refl_fr)

sim.run(mp.at_beginning(mp.output_epsilon),
mp.to_appended("ez", mp.at_every(0.1, mp.output_efield_z)),
until=sim_duration)

sim.save_flux('refl-flux', refl)
sim.display_fluxes(refl)
import meep as mp
import metals
from run_parameters import *

cell=mp.Vector3(x_size, y_size)

sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df),
 component=mp.Ez,
 center=mp.Vector3(0, source_y),
 size=mp.Vector3(x_size, 0))]

geometry=[mp.Block(mp.Vector3(x_size, metal_thickness, 1e20),
   center=mp.Vector3(0, metal_y),
   material=metals.silver)]

pml_layers=[mp.PML(pmld)]

sim=mp.Simulation(cell_size=cell,
  boundary_layers=pml_layers,
  geometry=geometry,
  sources=sources,
  resolution=resolution)

refl_fr = mp.FluxRegion(center=mp.Vector3(0, fr_y), size=mp.Vector3(x_size/2-pmld, 0))

#pt=mp.Vector3(0, 0)

refl = sim.add_flux(fcen, df, nfreq, refl_fr)

sim.load_minus_flux('refl-flux', refl)

sim.run(mp.at_beginning(mp.output_epsilon),
mp.to_appended("ez", mp.at_every(0.1, mp.output_efield_z)),
until=sim_duration)

sim.display_fluxes(refl)
import meep as mp

susc_aluminium = [
mp.LorentzianSusceptibility(frequency=8.06e-21, gamma=0.0379, sigma=1.173e42),
mp.LorentzianSusceptibility(frequency=0.13066, gamma=0.2686, sigma=1940.97)
]

susc_silver = [
mp.LorentzianSusceptibility(frequency=1e-20, gamma=0.038715, sigma=4.4625e41),
mp.LorentzianSusceptibility(frequency=0.65815, gamma=3.1343, sigma=7.9247),
mp.LorentzianSusceptibility(frequency=3.6142, gamma=0.36456, sigma=0.50133),
mp.LorentzianSusceptibility(frequency=6.6017, gamma=0.052426, sigma=0.013329),
mp.LorentzianSusceptibility(frequency=7.3259, gamma=0.7388, sigma=0.82655),
mp.LorentzianSusceptibility(frequency=16.365, gamma=1.9511, sigma=1.1133)
]

silver = mp.Medium(epsilon=1.0, E_susceptibilities=susc_silver)
x_size=.6
y_size=.5
pmld=.10
metal_thickness=.1
sim_duration=11

fcen = 3.333
df = 3.0

source_y=-y_size/2+1.5*pmld
metal_y=y_size/2-metal_thickness/2-1.5*pmld
fr_y=((y_size/2-1.5*pmld-metal_thickness)+(1.5*pmld-y_size/2))/2

nfreq=250

resolution=398

if __name__=="__main__":
print("Simulation space from Y = ", -y_size/2, " to Y = ", y_size/2)
print("Source position = ", source_y)
print("Metal position = ", metal_y, " from ", metal_y-metal_thickness/2, " to ", metal_y+metal_thickness/2)
print("Flux measured at = ", fr_y)
___
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss