From f805934d483149e5aee0367e11bef4575bc09625 Mon Sep 17 00:00:00 2001 From: smartalecH Date: Mon, 24 Feb 2020 20:57:14 -0500 Subject: [PATCH 1/8] init --- python/meep.i | 1 + python/source.py | 61 +++++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 61 insertions(+), 1 deletion(-) diff --git a/python/meep.i b/python/meep.i index 3c8498c1e..b9f6c2ca3 100644 --- a/python/meep.i +++ b/python/meep.i @@ -1555,6 +1555,7 @@ PyObject *_get_array_slice_dimensions(meep::fields *f, const meep::volume &where from .source import ( ContinuousSource, CustomSource, + FilteredSource, EigenModeSource, GaussianSource, Source, diff --git a/python/source.py b/python/source.py index 889fa5cfe..3482e4434 100644 --- a/python/source.py +++ b/python/source.py @@ -2,7 +2,8 @@ import meep as mp from meep.geom import Vector3, check_nonnegative - +from scipy import signal +import numpy as np def check_positive(prop, val): if val > 0: @@ -80,6 +81,64 @@ def __init__(self, src_func, start_time=-1.0e20, end_time=1.0e20, center_frequen self.swigobj = mp.custom_src_time(src_func, start_time, end_time, center_frequency) self.swigobj.is_integrated = self.is_integrated +class FilteredSource(CustomSource): + def __init__(self,center_frequency,frequencies,frequency_response,num_taps,dt,time_src): + self.center_frequency=center_frequency + self.frequencies=frequencies + self.frequency_response=frequency_response + self.num_taps=500 + self.dt=dt + self.time_src=time_src + self.current_time = None + + f = self.func() + + # initialize super + super(FilteredSource, self).__init__(src_func=f,center_frequency=self.center_frequency) + + # calculate equivalent sample rate + self.fs = 1/self.dt + + # estimate impulse response from frequency response + self.estimate_impulse_response() + + def filter(self,t): + # shift feedforward memory + np.roll(self.memory, -1) + if self.current_time is None or self.current_time != t: + self.current_time = t + self.memory[0] = self.time_src.swigobj.dipole(t) + # calculate filter response + self.current_y = np.dot(self.memory,self.taps) + return self.current_y + + def estimate_impulse_response(self): + # calculate band edges from target frequencies + w = self.frequencies/(self.fs/2) * np.pi + D = self.frequency_response + self.taps = self.lstsqrs(self.num_taps,w,D) + + # allocate filter memory taps + self.memory = np.zeros(self.taps.shape,dtype=np.complex128) + + def func(self): + def _f(t): + return self.filter(t) + return _f + + def lstsqrs(self,num_taps,freqs,h_desired): + n_freqs = freqs.size + vandermonde = np.zeros((n_freqs,num_taps),dtype=np.complex128) + for iom, om in enumerate(freqs): + for it in range(num_taps): + vandermonde[iom,it] = np.exp(-1j*it*om) + + + a = np.matmul(np.linalg.pinv(vandermonde), h_desired) + _, h_hat = signal.freqz(a,worN=freqs) + l2_error = np.sqrt(np.sum(np.abs(h_hat - h_desired)**2)) + + return a class EigenModeSource(Source): From 20ae8530e268f11fae3d0efc4f6764c6b4cdddae Mon Sep 17 00:00:00 2001 From: smartalecH Date: Mon, 24 Feb 2020 21:27:27 -0500 Subject: [PATCH 2/8] updates --- python/source.py | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/python/source.py b/python/source.py index 3482e4434..fd9b6770e 100644 --- a/python/source.py +++ b/python/source.py @@ -103,13 +103,11 @@ def __init__(self,center_frequency,frequencies,frequency_response,num_taps,dt,ti self.estimate_impulse_response() def filter(self,t): - # shift feedforward memory - np.roll(self.memory, -1) if self.current_time is None or self.current_time != t: - self.current_time = t - self.memory[0] = self.time_src.swigobj.dipole(t) - # calculate filter response - self.current_y = np.dot(self.memory,self.taps) + self.current_time = t # increase current time (we went through all the current components) + np.roll(self.memory, -1) # shift feedforward memory + self.memory[0] = self.time_src.swigobj.dipole(t) # update current memory slot + self.current_y = np.dot(self.memory,self.taps) # calculate filter response as inner product of taps and memory return self.current_y def estimate_impulse_response(self): @@ -133,11 +131,9 @@ def lstsqrs(self,num_taps,freqs,h_desired): for it in range(num_taps): vandermonde[iom,it] = np.exp(-1j*it*om) - a = np.matmul(np.linalg.pinv(vandermonde), h_desired) _, h_hat = signal.freqz(a,worN=freqs) - l2_error = np.sqrt(np.sum(np.abs(h_hat - h_desired)**2)) - + self.l2_error = np.sqrt(np.sum(np.abs(h_hat - h_desired)**2)) return a class EigenModeSource(Source): From 58ae86a71dd9d951617964f989542f8b475dad33 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Wed, 26 Feb 2020 16:01:17 -0500 Subject: [PATCH 3/8] some updates --- python/source.py | 78 +++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 71 insertions(+), 7 deletions(-) diff --git a/python/source.py b/python/source.py index fd9b6770e..31ed2657a 100644 --- a/python/source.py +++ b/python/source.py @@ -4,6 +4,7 @@ from meep.geom import Vector3, check_nonnegative from scipy import signal import numpy as np +from scipy.interpolate import interp1d def check_positive(prop, val): if val > 0: @@ -86,8 +87,8 @@ def __init__(self,center_frequency,frequencies,frequency_response,num_taps,dt,ti self.center_frequency=center_frequency self.frequencies=frequencies self.frequency_response=frequency_response - self.num_taps=500 - self.dt=dt + self.num_taps=num_taps + self.dt=dt/2 # the "effective" dt needs double resolution for staggered yee grid self.time_src=time_src self.current_time = None @@ -103,18 +104,26 @@ def __init__(self,center_frequency,frequencies,frequency_response,num_taps,dt,ti self.estimate_impulse_response() def filter(self,t): + idx = int(np.round(t/self.dt)) + if idx >= self.num_taps: + return 0 + else: + return self.taps[idx] + + + '''print(idx) if self.current_time is None or self.current_time != t: self.current_time = t # increase current time (we went through all the current components) np.roll(self.memory, -1) # shift feedforward memory self.memory[0] = self.time_src.swigobj.dipole(t) # update current memory slot self.current_y = np.dot(self.memory,self.taps) # calculate filter response as inner product of taps and memory - return self.current_y + return self.current_y''' def estimate_impulse_response(self): # calculate band edges from target frequencies w = self.frequencies/(self.fs/2) * np.pi D = self.frequency_response - self.taps = self.lstsqrs(self.num_taps,w,D) + self.taps = self.spline_fit(self.num_taps,w,D) # allocate filter memory taps self.memory = np.zeros(self.taps.shape,dtype=np.complex128) @@ -124,16 +133,71 @@ def _f(t): return self.filter(t) return _f + def spline_fit(self,num_taps,freqs,h_desired): + num_taps = 2000 + # fit real part + real_x = np.concatenate(([0],freqs,[np.pi])) + real_y = np.concatenate(([0],np.real(h_desired),[0])) + fr = interp1d(real_x, real_y, kind='cubic') + + # fit imaginary part + imag_x = np.concatenate(([0],freqs,[np.pi])) + imag_y = np.concatenate(([0],np.imag(h_desired),[0])) + fi = interp1d(imag_x, imag_y, kind='cubic') + + # formulate hermitian filter response with specified number of taps + freqs_filter = numpy.fft.fftfreq#np.linspace(-np.pi,np.pi,num_taps) + zero_freq = np.argmin(np.abs(freqs_filter))+1 + print(freqs_filter[zero_freq]) + + filter_pos = fr(freqs_filter[zero_freq:]) + 1j*fi(freqs_filter[zero_freq:]) + filter_neg = np.flipud(np.conjugate(filter_pos)) + + if num_taps %2 == 0: #even + filter_both = np.concatenate((filter_pos,filter_neg)) + else: + filter_both = np.concatenate((filter_pos,filter_neg[:-1])) + + from matplotlib import pyplot as plt + plt.figure() + #plt.plot(freqs_filter,np.abs(filter_both)) + plt.plot(freqs_filter[zero_freq:],fr(freqs_filter[zero_freq:])) + plt.plot(freqs,np.real(h_desired),'o') + plt.show() + + print(num_taps) + print(filter_both.size) + quit() + + # ifft to get sampled impulse response + + return def lstsqrs(self,num_taps,freqs,h_desired): n_freqs = freqs.size - vandermonde = np.zeros((n_freqs,num_taps),dtype=np.complex128) + vandermonde_left = np.zeros((n_freqs,num_taps),dtype=np.complex128) + vandermonde_right = np.zeros((n_freqs,num_taps),dtype=np.complex128) for iom, om in enumerate(freqs): for it in range(num_taps): - vandermonde[iom,it] = np.exp(-1j*it*om) + vandermonde_left[iom,it] = np.exp(-1j*it*om) + vandermonde_right[iom,it] = np.exp(1j*it*om) + vandermonde = np.vstack((vandermonde_left,vandermonde_right)) + h_desired_full = np.hstack((h_desired,np.conj(h_desired))) - a = np.matmul(np.linalg.pinv(vandermonde), h_desired) + a = np.matmul(np.linalg.pinv(vandermonde), h_desired_full) _, h_hat = signal.freqz(a,worN=freqs) + + import matplotlib.pyplot as plt + plt.figure() + plt.plot(freqs,np.abs(h_desired)) + plt.plot(freqs,np.abs(h_hat),'--') + plt.show() + quit() + self.l2_error = np.sqrt(np.sum(np.abs(h_hat - h_desired)**2)) + print(self.l2_error) + + # account for dtft scaling + a = a*self.dt*np.sqrt(2) return a class EigenModeSource(Source): From 815b6e0b5506038083b6beecc8191bb682ea6acb Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Tue, 3 Mar 2020 16:29:22 -0500 Subject: [PATCH 4/8] fixed scaling --- python/source.py | 178 ++++++++++++++++++++++------------------------- src/meep.hpp | 6 -- tests/temp.py | 116 ++++++++++++++++++++++++++++++ 3 files changed, 199 insertions(+), 101 deletions(-) create mode 100644 tests/temp.py diff --git a/python/source.py b/python/source.py index 31ed2657a..352ab2067 100644 --- a/python/source.py +++ b/python/source.py @@ -83,122 +83,110 @@ def __init__(self, src_func, start_time=-1.0e20, end_time=1.0e20, center_frequen self.swigobj.is_integrated = self.is_integrated class FilteredSource(CustomSource): - def __init__(self,center_frequency,frequencies,frequency_response,num_taps,dt,time_src): + def __init__(self,center_frequency,frequencies,frequency_response,time_src,min_err=1e-6): self.center_frequency=center_frequency self.frequencies=frequencies - self.frequency_response=frequency_response - self.num_taps=num_taps - self.dt=dt/2 # the "effective" dt needs double resolution for staggered yee grid + + signal_fourier_transform = np.array([time_src.fourier_transform(f) for f in self.frequencies]) + self.frequency_response= signal_fourier_transform * frequency_response self.time_src=time_src self.current_time = None + self.min_err = min_err f = self.func() # initialize super - super(FilteredSource, self).__init__(src_func=f,center_frequency=self.center_frequency) - - # calculate equivalent sample rate - self.fs = 1/self.dt + super(FilteredSource, self).__init__(src_func=f,center_frequency=self.center_frequency,is_integrated=time_src.is_integrated) # estimate impulse response from frequency response self.estimate_impulse_response() - - def filter(self,t): - idx = int(np.round(t/self.dt)) - if idx >= self.num_taps: - return 0 - else: - return self.taps[idx] - - '''print(idx) - if self.current_time is None or self.current_time != t: - self.current_time = t # increase current time (we went through all the current components) - np.roll(self.memory, -1) # shift feedforward memory - self.memory[0] = self.time_src.swigobj.dipole(t) # update current memory slot - self.current_y = np.dot(self.memory,self.taps) # calculate filter response as inner product of taps and memory - return self.current_y''' - - def estimate_impulse_response(self): - # calculate band edges from target frequencies - w = self.frequencies/(self.fs/2) * np.pi - D = self.frequency_response - self.taps = self.spline_fit(self.num_taps,w,D) - # allocate filter memory taps - self.memory = np.zeros(self.taps.shape,dtype=np.complex128) + def __call__(self,t): + # simple RBF with gaussian kernel reduces to inner product at time step + return np.dot(self.gauss_t(t,self.frequencies,self.gaus_widths),self.nodes) def func(self): def _f(t): - return self.filter(t) + return self(t) return _f - def spline_fit(self,num_taps,freqs,h_desired): - num_taps = 2000 - # fit real part - real_x = np.concatenate(([0],freqs,[np.pi])) - real_y = np.concatenate(([0],np.real(h_desired),[0])) - fr = interp1d(real_x, real_y, kind='cubic') - - # fit imaginary part - imag_x = np.concatenate(([0],freqs,[np.pi])) - imag_y = np.concatenate(([0],np.imag(h_desired),[0])) - fi = interp1d(imag_x, imag_y, kind='cubic') - - # formulate hermitian filter response with specified number of taps - freqs_filter = numpy.fft.fftfreq#np.linspace(-np.pi,np.pi,num_taps) - zero_freq = np.argmin(np.abs(freqs_filter))+1 - print(freqs_filter[zero_freq]) - - filter_pos = fr(freqs_filter[zero_freq:]) + 1j*fi(freqs_filter[zero_freq:]) - filter_neg = np.flipud(np.conjugate(filter_pos)) - - if num_taps %2 == 0: #even - filter_both = np.concatenate((filter_pos,filter_neg)) - else: - filter_both = np.concatenate((filter_pos,filter_neg[:-1])) - - from matplotlib import pyplot as plt - plt.figure() - #plt.plot(freqs_filter,np.abs(filter_both)) - plt.plot(freqs_filter[zero_freq:],fr(freqs_filter[zero_freq:])) - plt.plot(freqs,np.real(h_desired),'o') - plt.show() - - print(num_taps) - print(filter_both.size) - quit() - - # ifft to get sampled impulse response - - return - def lstsqrs(self,num_taps,freqs,h_desired): - n_freqs = freqs.size - vandermonde_left = np.zeros((n_freqs,num_taps),dtype=np.complex128) - vandermonde_right = np.zeros((n_freqs,num_taps),dtype=np.complex128) - for iom, om in enumerate(freqs): - for it in range(num_taps): - vandermonde_left[iom,it] = np.exp(-1j*it*om) - vandermonde_right[iom,it] = np.exp(1j*it*om) - vandermonde = np.vstack((vandermonde_left,vandermonde_right)) - h_desired_full = np.hstack((h_desired,np.conj(h_desired))) + def estimate_impulse_response(self): + ''' + find gaussian weighting coefficients. + + TODO use optimizer to find optimal gaussian widths + ''' + # Use vandermonde matrix to calculate weights of each gaussian. + # Each gaussian is centered at each frequency point + h = self.frequency_response + def rbf_l2(fwidth): + vandermonde = np.zeros((self.frequencies.size,self.frequencies.size),dtype=np.complex128) + for ri, rf in enumerate(self.frequencies): + for ci, cf in enumerate(self.frequencies): + vandermonde[ri,ci] = self.gauss_f(rf,cf,fwidth) + + nodes = np.matmul(np.linalg.pinv(vandermonde),h) + h_hat = np.matmul(vandermonde,nodes) + l2_err = np.sum(np.abs(h-h_hat)**2) + return nodes, l2_err - a = np.matmul(np.linalg.pinv(vandermonde), h_desired_full) - _, h_hat = signal.freqz(a,worN=freqs) - - import matplotlib.pyplot as plt - plt.figure() - plt.plot(freqs,np.abs(h_desired)) - plt.plot(freqs,np.abs(h_hat),'--') - plt.show() - quit() + df = self.frequencies[2] - self.frequencies[1] + err_high = True + fwidth = 1/self.time_src.width - self.l2_error = np.sqrt(np.sum(np.abs(h_hat - h_desired)**2)) - print(self.l2_error) + # Iterate through smaller and smaller widths until error is small enough or width is distance between frequency points + while err_high: + nodes, l2_err = rbf_l2(fwidth) + if l2_err < self.min_err or fwidth < df: + err_high = False + else: + fwidth = 0.5 * fwidth + self.gaus_widths = fwidth + self.nodes = nodes - # account for dtft scaling - a = a*self.dt*np.sqrt(2) - return a + from matplotlib import pyplot as plt + + temp = self.gauss_f(self.frequencies[:,np.newaxis],self.frequencies,fwidth) + i_hat = np.inner(self.nodes,temp) + + '''plt.figure() + plt.subplot(2,1,1) + plt.title('') + plt.semilogy(self.frequencies,np.abs(h),label='Desired response') + plt.semilogy(self.frequencies,np.abs(i_hat),'--',label='Fit response') + plt.legend() + plt.ylabel('Magnitude') + plt.grid(True) + + plt.subplot(2,1,2) + plt.plot(self.frequencies,np.unwrap(np.angle(self.frequency_response)),label='Desired response') + plt.plot(self.frequencies,np.unwrap(np.angle(i_hat)),'--',label='Fit response') + plt.legend() + plt.ylabel('Phase') + plt.xlabel('Frequency') + plt.grid(True) + + plt.savefig('fit_results.png') + plt.show()''' + + def gauss_t(self,t,f0,fwidth): + s = 5 + w = 1.0 / fwidth + t0 = w * s + tt = (t - t0) + amp = 1.0 / (-1j*2*np.pi*f0) + return amp * np.exp(-tt * tt / (2 * w * w))*np.exp(-1j*2 * np.pi * f0 * tt) #/ np.sqrt(2*np.pi) + + def gauss_f(self,f,f0,fwidth): + s = 5 + w = 1.0 / fwidth + t0 = w * s + omega = 2.0 * np.pi * f + omega0 = 2.0 * np.pi * f0 + delta = (omega - omega0) * w + amp = self.center_frequency/f0#/np.sqrt(omega0)#1j / (omega0) + return amp * w * np.exp(1j*omega*t0) * np.exp(-0.5 * delta * delta) class EigenModeSource(Source): diff --git a/src/meep.hpp b/src/meep.hpp index 4d95e8289..e07f1f152 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -928,12 +928,6 @@ class custom_src_time : public src_time { : func(func), data(data), freq(f), start_time(float(st)), end_time(float(et)) {} virtual ~custom_src_time() {} - virtual std::complex current(double time, double dt) const { - if (is_integrated) - return src_time::current(time, dt); - else - return dipole(time); - } virtual std::complex dipole(double time) const { float rtime = float(time); if (rtime >= start_time && rtime <= end_time) diff --git a/tests/temp.py b/tests/temp.py new file mode 100644 index 000000000..60139892e --- /dev/null +++ b/tests/temp.py @@ -0,0 +1,116 @@ +import meep as mp +import numpy as np +import matplotlib.pyplot as plt +from scipy import signal +import meep_adjoint as mpa + +# --------------------------- # +# Design filter +# --------------------------- # + +freq_min = 1/1.7 +freq_max = 1/1.4 +nfreq = 200 +freqs = np.linspace(freq_min,freq_max,nfreq) +resolution = 10 +run_time = 4000 + +# --------------------------- # +# Run normal simulation +# --------------------------- # + +cell = [16,8,0] +geometry = [mp.Block(mp.Vector3(mp.inf,1,mp.inf), + center=mp.Vector3(), + material=mp.Medium(epsilon=12))] +fcen = 1/1.55 +gauss_src = mp.GaussianSource(frequency=fcen,fwidth=0.1/1.55,is_integrated=False) +sources = [mp.EigenModeSource(gauss_src, + eig_band=1, + size=[0,8], + center=[-6,0])] +pml_layers = [mp.PML(1.0)] + +sim = mp.Simulation(cell_size=cell, + boundary_layers=pml_layers, + geometry=geometry, + sources=sources, + resolution=resolution) +mon = sim.add_dft_fields([mp.Ez],freq_min,freq_max,nfreq,center=[0,0],size=[0,1]) + +sim.run(until=run_time) + +fields = np.zeros((nfreq,),dtype=np.complex128) +# just store one spatial point for each freq +for f in range(nfreq): + fields[f] = sim.get_dft_array(mon,mp.Ez,f)[10] + +# build simple bandpass filter +dt = sim.fields.dt +fs = 1/dt +freqs_scipy = freqs * np.pi / (fs/2) +num_taps = 321 +taps = signal.firwin(num_taps,[freq_min/(fs/2), freq_max/(fs/2)], pass_zero='bandstop',window='boxcar') + +w,h = signal.freqz(taps,worN=freqs_scipy) + +'''plt.figure() +plt.plot(w,np.abs(h)) +plt.show() +quit()''' + +# frequency domain calculation +desired_fields = h * fields + +# --------------------------- # +# Run filtered simulation +# --------------------------- # + +filtered_src = mp.FilteredSource(fcen,freqs,h,gauss_src) + +sources = [mp.EigenModeSource(filtered_src, + eig_band=1, + size=[0,8], + center=[-6,0])] + +sim.reset_meep() +sim.change_sources(sources) + +mon = sim.add_dft_fields([mp.Ez],freq_min,freq_max,nfreq,center=[0,0],size=[0,1]) +sim.run(until=run_time) + +fields_filtered = np.zeros((nfreq,),dtype=np.complex128) +# just store one spatial point for each freq +for f in range(nfreq): + fields_filtered[f] = sim.get_dft_array(mon,mp.Ez,f)[10] + +# --------------------------- # +# Compare results +# --------------------------- # +print(np.abs(fields_filtered/desired_fields), np.angle(fields_filtered/desired_fields)) + +plt.figure() +plt.subplot(2,1,1) +plt.semilogy(freqs,np.abs(desired_fields),label='Frequency Domain') +plt.semilogy(freqs,np.abs(fields_filtered),'-.',label='Time Domain') +#plt.plot(freqs,np.abs(fields),'--',label='Gaussian Src') +#plt.plot(freqs,np.abs(h),'--',label='Window') +plt.grid() +plt.xlabel('Frequency') +plt.ylabel('Magnitude') +plt.legend() + +plt.subplot(2,1,2) +plt.plot(freqs,np.unwrap(np.angle(desired_fields)),label='Frequency Domain') +plt.plot(freqs,np.unwrap(np.angle(fields_filtered)),'-.',label='Time Domain') +#plt.plot(freqs,np.unwrap(np.angle(fields)),'--',label='Gaussian Src') +#plt.plot(freqs,np.unwrap(np.angle(h)),'--',label='Window') +plt.grid() +plt.xlabel('Frequency') +plt.ylabel('Angle') +plt.legend() + +plt.tight_layout() +plt.savefig('filtered.png') + +plt.show() \ No newline at end of file From eeb1796194bfd7e100847008d4c91f4bc6772689 Mon Sep 17 00:00:00 2001 From: smartalecH Date: Wed, 4 Mar 2020 20:40:05 -0500 Subject: [PATCH 5/8] use dtft instead of dft --- python/source.py | 120 ++++++++++------------------------------------- src/meep.hpp | 6 +++ 2 files changed, 30 insertions(+), 96 deletions(-) diff --git a/python/source.py b/python/source.py index 352ab2067..9aa1fd585 100644 --- a/python/source.py +++ b/python/source.py @@ -4,7 +4,7 @@ from meep.geom import Vector3, check_nonnegative from scipy import signal import numpy as np -from scipy.interpolate import interp1d +from scipy.interpolate import Rbf, PchipInterpolator def check_positive(prop, val): if val > 0: @@ -83,110 +83,38 @@ def __init__(self, src_func, start_time=-1.0e20, end_time=1.0e20, center_frequen self.swigobj.is_integrated = self.is_integrated class FilteredSource(CustomSource): - def __init__(self,center_frequency,frequencies,frequency_response,time_src,min_err=1e-6): + def __init__(self,center_frequency,frequencies,frequency_response,dt,T,time_src,min_err=1e-6): self.center_frequency=center_frequency self.frequencies=frequencies - - signal_fourier_transform = np.array([time_src.fourier_transform(f) for f in self.frequencies]) - self.frequency_response= signal_fourier_transform * frequency_response self.time_src=time_src - self.current_time = None self.min_err = min_err - f = self.func() + # calculate dtft of input signal + signal_t = np.array([time_src.swigobj.current(t,dt) for t in np.arange(0,T,dt)]) + signal_dtft = np.zeros((frequencies.shape),dtype=np.complex128) + for n, st in enumerate(signal_t): + for fi, f in enumerate(frequencies): + signal_dtft[fi] += np.exp(1j*2*np.pi*f*n*dt)*st - # initialize super - super(FilteredSource, self).__init__(src_func=f,center_frequency=self.center_frequency,is_integrated=time_src.is_integrated) + # multiply sampled dft of input signal with filter transfer function + H = signal_dtft #* frequency_response - # estimate impulse response from frequency response - self.estimate_impulse_response() - + # fit final frequency response to rbf + H_rbf = Rbf(frequencies, H, function='gaussian') + f_rbf = np.arange(0,1/dt,1/T) - def __call__(self,t): - # simple RBF with gaussian kernel reduces to inner product at time step - return np.dot(self.gauss_t(t,self.frequencies,self.gaus_widths),self.nodes) - - def func(self): - def _f(t): - return self(t) - return _f - - def estimate_impulse_response(self): - ''' - find gaussian weighting coefficients. - - TODO use optimizer to find optimal gaussian widths - ''' - # Use vandermonde matrix to calculate weights of each gaussian. - # Each gaussian is centered at each frequency point - h = self.frequency_response - def rbf_l2(fwidth): - vandermonde = np.zeros((self.frequencies.size,self.frequencies.size),dtype=np.complex128) - for ri, rf in enumerate(self.frequencies): - for ci, cf in enumerate(self.frequencies): - vandermonde[ri,ci] = self.gauss_f(rf,cf,fwidth) - - nodes = np.matmul(np.linalg.pinv(vandermonde),h) - h_hat = np.matmul(vandermonde,nodes) - l2_err = np.sum(np.abs(h-h_hat)**2) - return nodes, l2_err - - df = self.frequencies[2] - self.frequencies[1] - err_high = True - fwidth = 1/self.time_src.width + # estimate impulse response of final frequency response using ifft + h = np.flipud(np.fft.ifft(H_rbf(f_rbf))) # flip signal becuase fft convention is backwards of meeps + t_h = np.arange(0,T,dt) + + # fit impulse response to function to make implementation easy + h_rbf = PchipInterpolator(t_h, h) # we don't need a ton of extra accuracy at this point -- just speed + + def f_temp(t): + return np.asscalar(h_rbf(t)) - # Iterate through smaller and smaller widths until error is small enough or width is distance between frequency points - while err_high: - nodes, l2_err = rbf_l2(fwidth) - if l2_err < self.min_err or fwidth < df: - err_high = False - else: - fwidth = 0.5 * fwidth - self.gaus_widths = fwidth - self.nodes = nodes - - from matplotlib import pyplot as plt - - temp = self.gauss_f(self.frequencies[:,np.newaxis],self.frequencies,fwidth) - i_hat = np.inner(self.nodes,temp) - - '''plt.figure() - plt.subplot(2,1,1) - plt.title('') - plt.semilogy(self.frequencies,np.abs(h),label='Desired response') - plt.semilogy(self.frequencies,np.abs(i_hat),'--',label='Fit response') - plt.legend() - plt.ylabel('Magnitude') - plt.grid(True) - - plt.subplot(2,1,2) - plt.plot(self.frequencies,np.unwrap(np.angle(self.frequency_response)),label='Desired response') - plt.plot(self.frequencies,np.unwrap(np.angle(i_hat)),'--',label='Fit response') - plt.legend() - plt.ylabel('Phase') - plt.xlabel('Frequency') - plt.grid(True) - - plt.savefig('fit_results.png') - plt.show()''' - - def gauss_t(self,t,f0,fwidth): - s = 5 - w = 1.0 / fwidth - t0 = w * s - tt = (t - t0) - amp = 1.0 / (-1j*2*np.pi*f0) - return amp * np.exp(-tt * tt / (2 * w * w))*np.exp(-1j*2 * np.pi * f0 * tt) #/ np.sqrt(2*np.pi) - - def gauss_f(self,f,f0,fwidth): - s = 5 - w = 1.0 / fwidth - t0 = w * s - omega = 2.0 * np.pi * f - omega0 = 2.0 * np.pi * f0 - delta = (omega - omega0) * w - amp = self.center_frequency/f0#/np.sqrt(omega0)#1j / (omega0) - return amp * w * np.exp(1j*omega*t0) * np.exp(-0.5 * delta * delta) + # initialize super + super(FilteredSource, self).__init__(src_func=f_temp,center_frequency=self.center_frequency,is_integrated=time_src.is_integrated) class EigenModeSource(Source): diff --git a/src/meep.hpp b/src/meep.hpp index e07f1f152..4d95e8289 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -928,6 +928,12 @@ class custom_src_time : public src_time { : func(func), data(data), freq(f), start_time(float(st)), end_time(float(et)) {} virtual ~custom_src_time() {} + virtual std::complex current(double time, double dt) const { + if (is_integrated) + return src_time::current(time, dt); + else + return dipole(time); + } virtual std::complex dipole(double time) const { float rtime = float(time); if (rtime >= start_time && rtime <= end_time) From 82d056f3c1ef2e669311aba9af93282ddf21abd0 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Thu, 5 Mar 2020 13:23:43 -0500 Subject: [PATCH 6/8] do inverse dtft --- python/source.py | 98 ++++++++++++++++++------------------------------ src/meep.hpp | 6 +++ tests/temp.py | 14 +++---- 3 files changed, 49 insertions(+), 69 deletions(-) diff --git a/python/source.py b/python/source.py index 352ab2067..87f39430a 100644 --- a/python/source.py +++ b/python/source.py @@ -4,7 +4,7 @@ from meep.geom import Vector3, check_nonnegative from scipy import signal import numpy as np -from scipy.interpolate import interp1d +from scipy.special import erf def check_positive(prop, val): if val > 0: @@ -83,35 +83,51 @@ def __init__(self, src_func, start_time=-1.0e20, end_time=1.0e20, center_frequen self.swigobj.is_integrated = self.is_integrated class FilteredSource(CustomSource): - def __init__(self,center_frequency,frequencies,frequency_response,time_src,min_err=1e-6): + def __init__(self,center_frequency,frequencies,frequency_response,dt,T,time_src,min_err=1e-6): + dt = dt/2 + self.dt = dt self.center_frequency=center_frequency self.frequencies=frequencies - - signal_fourier_transform = np.array([time_src.fourier_transform(f) for f in self.frequencies]) - self.frequency_response= signal_fourier_transform * frequency_response self.time_src=time_src - self.current_time = None self.min_err = min_err - f = self.func() - # initialize super - super(FilteredSource, self).__init__(src_func=f,center_frequency=self.center_frequency,is_integrated=time_src.is_integrated) + # calculate dtft of input signal + signal_t = np.array([time_src.swigobj.current(t,dt) for t in np.arange(0,T,dt)]) # time domain signal + signal_dtft = np.exp(1j*2*np.pi*frequencies[:,np.newaxis]*np.arange(0,signal_t.size)[np.newaxis,:]*dt)@signal_t # vectorize dtft for speed - # estimate impulse response from frequency response - self.estimate_impulse_response() - + # multiply sampled dft of input signal with filter transfer function + H = signal_dtft#* frequency_response + + self.estimate_impulse_response(H) + # initialize super + super(FilteredSource, self).__init__(src_func=f,center_frequency=self.center_frequency,is_integrated=False) + + def gaussian(self,f,f0,fwidth): + return np.exp(-0.5*((f-f0)/fwidth)**2) + def antiderivative(self,f,n,f0,fwidth,T): + a = np.sqrt(np.pi/2)*fwidth + phase = np.exp(-2*np.pi*n*T*(np.pi*n*T*fwidth*fwidth + 1j*f0)) + kernel = erf(f/(np.sqrt(2)*fwidth)-f0/(np.sqrt(2)*fwidth)+1j*np.sqrt(2)*np.pi*n*fwidth*T) + return a*phase*kernel + def dtft_gaussian(self,n,f0,fwidth,T): + f_start = 0 + f_end = 1/T + return T*(self.antiderivative(f_end,n,f0,fwidth,T) - self.antiderivative(f_start,n,f0,fwidth,T)) def __call__(self,t): + n = int(np.round(t/self.dt)) + #print(t/self.dt,n) + vec = self.dtft_gaussian(n,self.frequencies,self.gaus_widths,self.dt) # simple RBF with gaussian kernel reduces to inner product at time step - return np.dot(self.gauss_t(t,self.frequencies,self.gaus_widths),self.nodes) + return np.dot(vec,self.nodes) def func(self): def _f(t): return self(t) return _f - def estimate_impulse_response(self): + def estimate_impulse_response(self,H): ''' find gaussian weighting coefficients. @@ -119,16 +135,11 @@ def estimate_impulse_response(self): ''' # Use vandermonde matrix to calculate weights of each gaussian. # Each gaussian is centered at each frequency point - h = self.frequency_response def rbf_l2(fwidth): - vandermonde = np.zeros((self.frequencies.size,self.frequencies.size),dtype=np.complex128) - for ri, rf in enumerate(self.frequencies): - for ci, cf in enumerate(self.frequencies): - vandermonde[ri,ci] = self.gauss_f(rf,cf,fwidth) - - nodes = np.matmul(np.linalg.pinv(vandermonde),h) - h_hat = np.matmul(vandermonde,nodes) - l2_err = np.sum(np.abs(h-h_hat)**2) + vandermonde = self.gaussian(self.frequencies[:,np.newaxis],self.frequencies[np.newaxis,:],fwidth) + nodes = np.matmul(np.linalg.pinv(vandermonde),H) + H_hat = np.matmul(vandermonde,nodes) + l2_err = np.sum(np.abs(H-H_hat)**2) return nodes, l2_err df = self.frequencies[2] - self.frequencies[1] @@ -142,52 +153,15 @@ def rbf_l2(fwidth): err_high = False else: fwidth = 0.5 * fwidth + print(l2_err) self.gaus_widths = fwidth self.nodes = nodes from matplotlib import pyplot as plt - temp = self.gauss_f(self.frequencies[:,np.newaxis],self.frequencies,fwidth) + temp = self.gaussian(self.frequencies[:,np.newaxis],self.frequencies,fwidth) i_hat = np.inner(self.nodes,temp) - '''plt.figure() - plt.subplot(2,1,1) - plt.title('') - plt.semilogy(self.frequencies,np.abs(h),label='Desired response') - plt.semilogy(self.frequencies,np.abs(i_hat),'--',label='Fit response') - plt.legend() - plt.ylabel('Magnitude') - plt.grid(True) - - plt.subplot(2,1,2) - plt.plot(self.frequencies,np.unwrap(np.angle(self.frequency_response)),label='Desired response') - plt.plot(self.frequencies,np.unwrap(np.angle(i_hat)),'--',label='Fit response') - plt.legend() - plt.ylabel('Phase') - plt.xlabel('Frequency') - plt.grid(True) - - plt.savefig('fit_results.png') - plt.show()''' - - def gauss_t(self,t,f0,fwidth): - s = 5 - w = 1.0 / fwidth - t0 = w * s - tt = (t - t0) - amp = 1.0 / (-1j*2*np.pi*f0) - return amp * np.exp(-tt * tt / (2 * w * w))*np.exp(-1j*2 * np.pi * f0 * tt) #/ np.sqrt(2*np.pi) - - def gauss_f(self,f,f0,fwidth): - s = 5 - w = 1.0 / fwidth - t0 = w * s - omega = 2.0 * np.pi * f - omega0 = 2.0 * np.pi * f0 - delta = (omega - omega0) * w - amp = self.center_frequency/f0#/np.sqrt(omega0)#1j / (omega0) - return amp * w * np.exp(1j*omega*t0) * np.exp(-0.5 * delta * delta) - class EigenModeSource(Source): def __init__(self, diff --git a/src/meep.hpp b/src/meep.hpp index e07f1f152..4d95e8289 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -928,6 +928,12 @@ class custom_src_time : public src_time { : func(func), data(data), freq(f), start_time(float(st)), end_time(float(et)) {} virtual ~custom_src_time() {} + virtual std::complex current(double time, double dt) const { + if (is_integrated) + return src_time::current(time, dt); + else + return dipole(time); + } virtual std::complex dipole(double time) const { float rtime = float(time); if (rtime >= start_time && rtime <= end_time) diff --git a/tests/temp.py b/tests/temp.py index 60139892e..33d1483f8 100644 --- a/tests/temp.py +++ b/tests/temp.py @@ -10,10 +10,10 @@ freq_min = 1/1.7 freq_max = 1/1.4 -nfreq = 200 +nfreq = 300 freqs = np.linspace(freq_min,freq_max,nfreq) resolution = 10 -run_time = 4000 +run_time = 400 # --------------------------- # # Run normal simulation @@ -66,7 +66,7 @@ # Run filtered simulation # --------------------------- # -filtered_src = mp.FilteredSource(fcen,freqs,h,gauss_src) +filtered_src = mp.FilteredSource(fcen,freqs,h,dt,run_time,gauss_src) sources = [mp.EigenModeSource(filtered_src, eig_band=1, @@ -87,11 +87,11 @@ # --------------------------- # # Compare results # --------------------------- # -print(np.abs(fields_filtered/desired_fields), np.angle(fields_filtered/desired_fields)) +print(np.abs(fields_filtered/fields), np.angle(fields_filtered/fields)) plt.figure() plt.subplot(2,1,1) -plt.semilogy(freqs,np.abs(desired_fields),label='Frequency Domain') +plt.semilogy(freqs,np.abs(fields),label='Frequency Domain') plt.semilogy(freqs,np.abs(fields_filtered),'-.',label='Time Domain') #plt.plot(freqs,np.abs(fields),'--',label='Gaussian Src') #plt.plot(freqs,np.abs(h),'--',label='Window') @@ -101,7 +101,7 @@ plt.legend() plt.subplot(2,1,2) -plt.plot(freqs,np.unwrap(np.angle(desired_fields)),label='Frequency Domain') +plt.plot(freqs,np.unwrap(np.angle(fields)),label='Frequency Domain') plt.plot(freqs,np.unwrap(np.angle(fields_filtered)),'-.',label='Time Domain') #plt.plot(freqs,np.unwrap(np.angle(fields)),'--',label='Gaussian Src') #plt.plot(freqs,np.unwrap(np.angle(h)),'--',label='Window') @@ -111,6 +111,6 @@ plt.legend() plt.tight_layout() -plt.savefig('filtered.png') +plt.savefig('filtered_{}.png'.format(resolution)) plt.show() \ No newline at end of file From c9775396af7e043ff1e34855fce86bfa72a2f305 Mon Sep 17 00:00:00 2001 From: smartalecH Date: Wed, 11 Mar 2020 21:14:21 -0400 Subject: [PATCH 7/8] use sinc instead of gaussian --- python/source.py | 67 +++++++++++++++++++++++++++++++----------------- 1 file changed, 43 insertions(+), 24 deletions(-) diff --git a/python/source.py b/python/source.py index 9aa1fd585..d6ed33888 100644 --- a/python/source.py +++ b/python/source.py @@ -67,6 +67,7 @@ def __init__(self, frequency=None, width=0, fwidth=float('inf'), start_time=0, c self.swigobj = mp.gaussian_src_time(self.frequency, self.width, self.start_time, self.start_time + 2 * self.width * self.cutoff) self.swigobj.is_integrated = self.is_integrated + self.peak_time = 0.5 * (self.start_time + self.start_time + 2 * self.width * self.cutoff) def fourier_transform(self, freq): return self.swigobj.fourier_transform(freq) @@ -83,38 +84,56 @@ def __init__(self, src_func, start_time=-1.0e20, end_time=1.0e20, center_frequen self.swigobj.is_integrated = self.is_integrated class FilteredSource(CustomSource): - def __init__(self,center_frequency,frequencies,frequency_response,dt,T,time_src,min_err=1e-6): - self.center_frequency=center_frequency + def __init__(self,center_frequency,frequencies,frequency_response,dt,T,time_src): + dt = dt/2 # divide by two to compensate for staggered E,H time interval + self.dt = dt self.frequencies=frequencies - self.time_src=time_src - self.min_err = min_err + self.N = np.round(T/dt) + f = self.func() # calculate dtft of input signal - signal_t = np.array([time_src.swigobj.current(t,dt) for t in np.arange(0,T,dt)]) - signal_dtft = np.zeros((frequencies.shape),dtype=np.complex128) - for n, st in enumerate(signal_t): - for fi, f in enumerate(frequencies): - signal_dtft[fi] += np.exp(1j*2*np.pi*f*n*dt)*st + signal_t = np.array([time_src.swigobj.current(t,dt) for t in np.arange(0,T,dt)]) # time domain signal + signal_dtft = np.exp(1j*2*np.pi*frequencies[:,np.newaxis]*np.arange(0,signal_t.size)[np.newaxis,:]*dt)@signal_t # vectorize dtft for speed # multiply sampled dft of input signal with filter transfer function - H = signal_dtft #* frequency_response + H = signal_dtft * frequency_response - # fit final frequency response to rbf - H_rbf = Rbf(frequencies, H, function='gaussian') - f_rbf = np.arange(0,1/dt,1/T) + # estimate the impulse response using a sinc function RBN + self.nodes, self.err = self.estimate_impulse_response(H) - # estimate impulse response of final frequency response using ifft - h = np.flipud(np.fft.ifft(H_rbf(f_rbf))) # flip signal becuase fft convention is backwards of meeps - t_h = np.arange(0,T,dt) - - # fit impulse response to function to make implementation easy - h_rbf = PchipInterpolator(t_h, h) # we don't need a ton of extra accuracy at this point -- just speed - - def f_temp(t): - return np.asscalar(h_rbf(t)) - # initialize super - super(FilteredSource, self).__init__(src_func=f_temp,center_frequency=self.center_frequency,is_integrated=time_src.is_integrated) + super(FilteredSource, self).__init__(src_func=f,center_frequency=center_frequency,is_integrated=False) + + def sinc(self,f,f0): + omega = 2*np.pi*f + omega0 = 2*np.pi*f0 + num = np.where(f == f0, self.N, (1-np.exp(1j*(self.N+1)*(omega-omega0)*self.dt))) + den = np.where(f == f0, 1, (1-np.exp(1j*(omega-omega0)*self.dt))) + return num/den + + def rect(self,n,f0): + omega0 = 2*np.pi*f0 + #return np.exp(-1j*omega0*(n-self.N/2)*self.dt) + return np.where(n < 0 or n > self.N,0,np.exp(-1j*omega0*(n)*self.dt)) + + def __call__(self,t): + n = int(np.round((t)/self.dt)) + vec = self.rect(n,self.frequencies) + return np.dot(vec,self.nodes) + + def func(self): + def _f(t): + return self(t) + return _f + + def estimate_impulse_response(self,H): + # Use vandermonde matrix to calculate weights of each gaussian. + # Each sinc is centered at each frequency point + vandermonde = self.sinc(self.frequencies[:,np.newaxis],self.frequencies[np.newaxis,:]) + nodes = np.matmul(np.linalg.pinv(vandermonde),H) + H_hat = np.matmul(vandermonde,nodes) + l2_err = np.sum(np.abs(H-H_hat)**2/np.abs(H)**2) + return nodes, l2_err class EigenModeSource(Source): From 4f382ab8e09e1bace34853bec382ddfa4b123def Mon Sep 17 00:00:00 2001 From: smartalecH Date: Fri, 20 Mar 2020 14:46:46 -0400 Subject: [PATCH 8/8] success --- python/source.py | 49 ++++++++++++++++++++++++++++++++++++++++-------- tests/temp.py | 26 ++++++++----------------- 2 files changed, 49 insertions(+), 26 deletions(-) diff --git a/python/source.py b/python/source.py index 7f7ed1bf7..7dc941726 100644 --- a/python/source.py +++ b/python/source.py @@ -103,23 +103,56 @@ def __init__(self,center_frequency,frequencies,frequency_response,dt,T,time_src) self.nodes, self.err = self.estimate_impulse_response(H) # initialize super - super(FilteredSource, self).__init__(src_func=f,center_frequency=center_frequency,is_integrated=False) + super(FilteredSource, self).__init__(end_time=T,src_func=f,center_frequency=center_frequency,is_integrated=False) + def cos_window_td(self,a,n,f0): + omega0 = 2*np.pi*f0 + cos_sum = np.sum([(-1)**k * a[k] * np.cos(2*np.pi*n*k/self.N )for k in range(len(a))]) + return np.where(n < 0 or n > self.N,0, + np.exp(-1j*omega0*n*self.dt) * (cos_sum) + ) + + def cos_window_fd(self,a,f,f0): + omega = 2*np.pi*f + omega0 = 2*np.pi*f0 + df = 1/(self.N*self.dt) + cos_sum = np.zeros((f*f0).shape,dtype=np.complex128)#a[0]*self.sinc(f,f0) + for k in range(len(a)): + cos_sum += (-1)**k * a[k]/2 * self.sinc(f,f0-k*df) + (-1)**k * a[k]/2 * self.sinc(f,f0+k*df) + return cos_sum + def sinc(self,f,f0): omega = 2*np.pi*f omega0 = 2*np.pi*f0 - num = np.where(f == f0, self.N, (1-np.exp(1j*(self.N+1)*(omega-omega0)*self.dt))) - den = np.where(f == f0, 1, (1-np.exp(1j*(omega-omega0)*self.dt))) + num = np.where(omega == omega0, self.N, (1-np.exp(1j*(self.N+1)*(omega-omega0)*self.dt))) + den = np.where(omega == omega0, 1, (1-np.exp(1j*(omega-omega0)*self.dt))) return num/den - def rect(self,n,f0): + def rect(self,t,f0): + n = int(np.round((t)/self.dt)) omega0 = 2*np.pi*f0 - #return np.exp(-1j*omega0*(n-self.N/2)*self.dt) return np.where(n < 0 or n > self.N,0,np.exp(-1j*omega0*(n)*self.dt)) + + def hann(self,t,f0): + n = int(np.round((t)/self.dt)) + a = [0.5, 0.5] + return self.cos_window_td(a,n,f0) + + def hann_dtft(self,f,f0): + a = [0.5, 0.5] + return self.cos_window_fd(a,f,f0) + + def nuttall(self,t,f0): + n = int(np.round((t)/self.dt)) + a = [0.355768, 0.4873960, 0.144232, 0.012604] + return self.cos_window_td(a,n,f0) + + def nutall_dtft(self,f,f0): + a = [0.355768, 0.4873960, 0.144232, 0.012604] + return self.cos_window_fd(a,f,f0) def __call__(self,t): - n = int(np.round((t)/self.dt)) - vec = self.rect(n,self.frequencies) + vec = self.nuttall(t,self.frequencies) return np.dot(vec,self.nodes) def func(self): @@ -130,7 +163,7 @@ def _f(t): def estimate_impulse_response(self,H): # Use vandermonde matrix to calculate weights of each gaussian. # Each sinc is centered at each frequency point - vandermonde = self.sinc(self.frequencies[:,np.newaxis],self.frequencies[np.newaxis,:]) + vandermonde = self.nutall_dtft(self.frequencies[:,np.newaxis],self.frequencies[np.newaxis,:]) nodes = np.matmul(np.linalg.pinv(vandermonde),H) H_hat = np.matmul(vandermonde,nodes) l2_err = np.sum(np.abs(H-H_hat)**2/np.abs(H)**2) diff --git a/tests/temp.py b/tests/temp.py index 33d1483f8..8bb314076 100644 --- a/tests/temp.py +++ b/tests/temp.py @@ -10,10 +10,10 @@ freq_min = 1/1.7 freq_max = 1/1.4 -nfreq = 300 +nfreq = 100 freqs = np.linspace(freq_min,freq_max,nfreq) resolution = 10 -run_time = 400 +run_time = 500 # --------------------------- # # Run normal simulation @@ -38,7 +38,7 @@ resolution=resolution) mon = sim.add_dft_fields([mp.Ez],freq_min,freq_max,nfreq,center=[0,0],size=[0,1]) -sim.run(until=run_time) +sim.run(until_after_sources=run_time) fields = np.zeros((nfreq,),dtype=np.complex128) # just store one spatial point for each freq @@ -54,19 +54,14 @@ w,h = signal.freqz(taps,worN=freqs_scipy) -'''plt.figure() -plt.plot(w,np.abs(h)) -plt.show() -quit()''' - # frequency domain calculation desired_fields = h * fields # --------------------------- # # Run filtered simulation # --------------------------- # - -filtered_src = mp.FilteredSource(fcen,freqs,h,dt,run_time,gauss_src) +T = 1200 +filtered_src = mp.FilteredSource(fcen,freqs,h,dt,T,gauss_src) sources = [mp.EigenModeSource(filtered_src, eig_band=1, @@ -77,7 +72,7 @@ sim.change_sources(sources) mon = sim.add_dft_fields([mp.Ez],freq_min,freq_max,nfreq,center=[0,0],size=[0,1]) -sim.run(until=run_time) +sim.run(until_after_sources=run_time) fields_filtered = np.zeros((nfreq,),dtype=np.complex128) # just store one spatial point for each freq @@ -87,24 +82,19 @@ # --------------------------- # # Compare results # --------------------------- # -print(np.abs(fields_filtered/fields), np.angle(fields_filtered/fields)) plt.figure() plt.subplot(2,1,1) -plt.semilogy(freqs,np.abs(fields),label='Frequency Domain') +plt.semilogy(freqs,np.abs(desired_fields),label='Frequency Domain') plt.semilogy(freqs,np.abs(fields_filtered),'-.',label='Time Domain') -#plt.plot(freqs,np.abs(fields),'--',label='Gaussian Src') -#plt.plot(freqs,np.abs(h),'--',label='Window') plt.grid() plt.xlabel('Frequency') plt.ylabel('Magnitude') plt.legend() plt.subplot(2,1,2) -plt.plot(freqs,np.unwrap(np.angle(fields)),label='Frequency Domain') +plt.plot(freqs,np.unwrap(np.angle(desired_fields)),label='Frequency Domain') plt.plot(freqs,np.unwrap(np.angle(fields_filtered)),'-.',label='Time Domain') -#plt.plot(freqs,np.unwrap(np.angle(fields)),'--',label='Gaussian Src') -#plt.plot(freqs,np.unwrap(np.angle(h)),'--',label='Window') plt.grid() plt.xlabel('Frequency') plt.ylabel('Angle')