diff --git a/src/dft.cpp b/src/dft.cpp index af5452fd1..1c5560cf0 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -33,8 +33,10 @@ namespace meep { struct dft_chunk_data { // for passing to field::loop_in_chunks as void* component c; int vc; - double omega_min, domega; + // TODO: change dft_chunk_data properties int Nomega; + std::vector omegas; + omegas.reserve(Nomega); complex stored_weight, extra_weight; double dt_factor; bool include_dV_and_interp_weights; @@ -84,9 +86,9 @@ dft_chunk::dft_chunk(fields_chunk *fc_, ivec is_, ivec ie_, vec s0_, vec s1_, ve sn = sn_; vc = data->vc; - omega_min = data->omega_min; - domega = data->domega; + // TODO: change dft_chunk properties Nomega = data->Nomega; + omegas.assign(data->omegas.begin(), data->omegas.end()); dft_phase = new complex[Nomega]; N = 1; @@ -105,6 +107,8 @@ dft_chunk::dft_chunk(fields_chunk *fc_, ivec is_, ivec ie_, vec s0_, vec s1_, ve dft_chunk::~dft_chunk() { delete[] dft; delete[] dft_phase; + // TODO: add omegas to destructor? + //delete[] omegas; // delete from fields_chunk list dft_chunk *cur = fc->dft_chunks; @@ -118,6 +122,8 @@ dft_chunk::~dft_chunk() { } void dft_flux::remove() { + // TODO: add freqs to remove? + // freqs.clear(); while (E) { dft_chunk *nxt = E->next_in_dft; delete E; @@ -145,11 +151,11 @@ static void add_dft_chunkloop(fields_chunk *fc, int ichunk, component cgrid, ive shift_phase * S.phase_shift(c, sn), shift, S, sn, chunkloop_data); } -dft_chunk *fields::add_dft(component c, const volume &where, double freq_min, double freq_max, - int Nfreq, bool include_dV_and_interp_weights, - complex stored_weight, dft_chunk *chunk_next, - bool sqrt_dV_and_interp_weights, complex extra_weight, - bool use_centered_grid, int vc) { +// TODO: overload fields::add_dft +dft_chunk *fields::add_dft(component c, const volume &where, std::vector freqs, int Nfreq, + bool include_dV_and_interp_weights, complex stored_weight, + dft_chunk *chunk_next, bool sqrt_dV_and_interp_weights, + complex extra_weight, bool use_centered_grid, int vc) { if (coordinate_mismatch(gv.dim, c)) return NULL; /* If you call add_dft before adding sources, it will do nothing @@ -164,10 +170,12 @@ dft_chunk *fields::add_dft(component c, const volume &where, double freq_min, do dft_chunk_data data; data.c = c; data.vc = vc; - if (Nfreq <= 1) freq_min = freq_max = (freq_min + freq_max) * 0.5; - data.omega_min = freq_min * 2 * pi; - data.domega = Nfreq <= 1 ? 0.0 : (freq_max * 2 * pi - data.omega_min) / (Nfreq - 1); + // TODO: change how freqs becomes omegas + if (Nfreq < 1) abort("Nfreq must be at least 1"); data.Nomega = Nfreq; + std::vector omegas; + for (int i = 0; i < Nfreq; i++) { omegas[i] = freqs[i] * 2 * pi; } + data.omegas.assign(omegas.begin(), omegas.end()); data.stored_weight = stored_weight; data.extra_weight = extra_weight; data.dt_factor = dt / sqrt(2.0 * pi); @@ -182,6 +190,36 @@ dft_chunk *fields::add_dft(component c, const volume &where, double freq_min, do return data.dft_chunks; } +dft_chunk *fields::add_dft(component c, const volume &where, double freq_min, double freq_max, + int Nfreq, bool include_dV_and_interp_weights, + complex stored_weight, dft_chunk *chunk_next, + bool sqrt_dV_and_interp_weights, complex extra_weight, + bool use_centered_grid, int vc) { + if (Nfreq < 1) abort("Nfreq must be at least 1"); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); + std::vector freqs; + for (int i = 0; i < Nfreq; i++) { + freqs[i] = freq_min + dfreq * i; + } + return add_dft(c, where, freqs, Nfreq, include_dV_and_interp_weights, stored_weight, chunk_next, + sqrt_dV_and_interp_weights, extra_weight, use_centered_grid, vc); +} + +// TODO: overload add_dft +dft_chunk *fields::add_dft(const volume_list *where, std::vector freqs, int Nfreq, + bool include_dV_and_interp_weights) { + dft_chunk *chunks = 0; + while (where) { + if (is_derived(where->c)) abort("derived_component invalid for dft"); + cdouble stored_weight = where->weight; + // TODO: change here + chunks = add_dft(component(where->c), where->v, freqs, Nfreq, include_dV_and_interp_weights, + stored_weight, chunks); + where = where->next; + } + return chunks; +} + dft_chunk *fields::add_dft(const volume_list *where, double freq_min, double freq_max, int Nfreq, bool include_dV_and_interp_weights) { dft_chunk *chunks = 0; @@ -195,6 +233,12 @@ dft_chunk *fields::add_dft(const volume_list *where, double freq_min, double fre return chunks; } +// TODO: overload add_dft_pt +dft_chunk *fields::add_dft_pt(component c, const vec &where, std::vector freqs, int Nfreq) { + // TODO: change here + return add_dft(c, where, freqs, Nfreq, false); +} + dft_chunk *fields::add_dft_pt(component c, const vec &where, double freq_min, double freq_max, int Nfreq) { return add_dft(c, where, freq_min, freq_max, Nfreq, false); @@ -218,7 +262,8 @@ void dft_chunk::update_dft(double time) { if (!fc->f[c][0]) return; for (int i = 0; i < Nomega; ++i) - dft_phase[i] = polar(1.0, (omega_min + i * domega) * time) * scale; + // TODO: change here, check this whole function + dft_phase[i] = polar(1.0, omegas[i] * time) * scale; int numcmp = fc->f[c][1] ? 2 : 1; @@ -246,11 +291,13 @@ void dft_chunk::update_dft(double time) { if (numcmp == 2) { complex fc(f[0], f[1]); for (int i = 0; i < Nomega; ++i) + // TODO: check this dft[Nomega * idx_dft + i] += dft_phase[i] * fc; } else { realnum fr = f[0]; for (int i = 0; i < Nomega; ++i) + // TODO: check this dft[Nomega * idx_dft + i] += dft_phase[i] * fr; } idx_dft++; @@ -335,20 +382,35 @@ void load_dft_hdf5(dft_chunk *dft_chunks, component c, h5file *file, const char load_dft_hdf5(dft_chunks, component_name(c), file, dprefix); } +// TODO: overload dft_flux constructor +dft_flux::dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, + std::vector fs, int Nf, const volume &where_, direction normal_direction_, + bool use_symmetry_) + : Nfreq(Nf), E(E_), H(H_), cE(cE_), cH(cH_), where(where_), normal_direction(normal_direction_), + use_symmetry(use_symmetry_) { + // TODO: change here + if (Nf < 1) abort("Nf must be at least 1"); + freqs.assign(fs.begin(), fs.end()); +} + dft_flux::dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, double fmin, double fmax, int Nf, const volume &where_, direction normal_direction_, bool use_symmetry_) : Nfreq(Nf), E(E_), H(H_), cE(cE_), cH(cH_), where(where_), normal_direction(normal_direction_), use_symmetry(use_symmetry_) { - if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5; - freq_min = fmin; - dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); + if (Nf < 1) abort("Nf must be at least 1"); + double dfreqs = (fmax - fmin) / (Nf - 1); + std::vector fs; + for (int i = 0; i < Nf; i++) { + fs[i] = fmin + dfreqs * i; + } + freqs.assign(fs.begin(), fs.end()); } dft_flux::dft_flux(const dft_flux &f) : where(f.where) { - freq_min = f.freq_min; + // TODO: change here Nfreq = f.Nfreq; - dfreq = f.dfreq; + freqs.assign(f.freqs.begin(), f.freqs.end()); E = f.E; H = f.H; cE = f.cE; @@ -401,10 +463,12 @@ void dft_flux::scale_dfts(complex scale) { if (H) H->scale_dft(scale); } -dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double freq_max, - int Nfreq, bool use_symmetry) { +// TODO: overload fields::add_dft_flux +dft_flux fields::add_dft_flux(const volume_list *where_, std::vector freqs, int Nfreq, + bool use_symmetry) { if (!where_) // handle empty list of volumes - return dft_flux(Ex, Hy, NULL, NULL, freq_min, freq_max, Nfreq, v, NO_DIRECTION, use_symmetry); + // TODO: change here + return dft_flux(Ex, Hy, NULL, NULL, freqs, Nfreq, v, NO_DIRECTION, use_symmetry); dft_chunk *E = 0, *H = 0; component cE[2] = {Ex, Ey}, cH[2] = {Hy, Hx}; @@ -437,9 +501,10 @@ dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double } for (int i = 0; i < 2; ++i) { - E = add_dft(cE[i], where->v, freq_min, freq_max, Nfreq, true, - where->weight * double(1 - 2 * i), E); - H = add_dft(cH[i], where->v, freq_min, freq_max, Nfreq, false, 1.0, H); + // TODO: change here + E = add_dft(cE[i], where->v, freqs, Nfreq, true, where->weight * double(1 - 2 * i), E); + // TODO: change here + H = add_dft(cH[i], where->v, freqs, Nfreq, false, 1.0, H); } where = where->next; @@ -449,21 +514,47 @@ dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double // if the volume list has only one entry, store its component's direction. // if the volume list has > 1 entry, store NO_DIRECTION. direction flux_dir = (where_->next ? NO_DIRECTION : component_direction(where_->c)); - return dft_flux(cE[0], cH[0], E, H, freq_min, freq_max, Nfreq, firstvol, flux_dir, use_symmetry); + // TODO: change here + return dft_flux(cE[0], cH[0], E, H, freqs, Nfreq, firstvol, flux_dir, use_symmetry); +} + +dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double freq_max, + int Nfreq, bool use_symmetry) { + if (Nfreq < 1) abort("Nfreq must be at least 1"); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); + std::vector freqs; + for (int i = 0; i < Nfreq; i++) { + freqs[i] = freq_min + dfreq * i; + } + return add_dft_flux(where_, freqs, Nfreq, use_symmetry); +} + +// TODO: overload dft_energy constructor +dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, + std::vector fs, int Nf, const volume &where_) + : Nfreq(Nf), E(E_), H(H_), D(D_), B(B_), where(where_) { + // TODO: change here + if (Nf < 1) abort("Nf must be at least 1"); + freqs.assign(fs.begin(), fs.end()); } dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double fmin, double fmax, int Nf, const volume &where_) : Nfreq(Nf), E(E_), H(H_), D(D_), B(B_), where(where_) { - if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5; - freq_min = fmin; - dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); + if (Nf < 1) abort("Nf must be at least 1"); + double dfreq = (fmax - fmin) / (Nf - 1); + std::vector fs; + for (int i = 0; i < Nf; i++) { + fs[i] = fmin + dfreq * i; + } + freqs.assign(fs.begin(), fs.end()); } +// TODO: change here dft_energy::dft_energy(const dft_energy &f) : where(f.where) { - freq_min = f.freq_min; + // TODO: change here Nfreq = f.Nfreq; - dfreq = f.dfreq; + freqs.assign(f.freqs.begin(), f.freqs.end()); E = f.E; H = f.H; D = f.D; @@ -511,11 +602,12 @@ double *dft_energy::total() { return F; } -dft_energy fields::add_dft_energy(const volume_list *where_, double freq_min, double freq_max, - int Nfreq) { +// TODO: overload fields::add_dft_energy +dft_energy fields::add_dft_energy(const volume_list *where_, std::vector freqs, int Nfreq) { if (!where_) // handle empty list of volumes - return dft_energy(NULL, NULL, NULL, NULL, freq_min, freq_max, Nfreq, v); + // TODO: change here + return dft_energy(NULL, NULL, NULL, NULL, freqs, Nfreq, v); dft_chunk *E = 0, *D = 0, *H = 0, *B = 0; volume firstvol(where_->v); @@ -523,16 +615,32 @@ dft_energy fields::add_dft_energy(const volume_list *where_, double freq_min, do volume_list *where_save = where; while (where) { LOOP_OVER_FIELD_DIRECTIONS(gv.dim, d) { - E = add_dft(direction_component(Ex, d), where->v, freq_min, freq_max, Nfreq, true, 1.0, E); - D = add_dft(direction_component(Dx, d), where->v, freq_min, freq_max, Nfreq, false, 1.0, D); - H = add_dft(direction_component(Hx, d), where->v, freq_min, freq_max, Nfreq, true, 1.0, H); - B = add_dft(direction_component(Bx, d), where->v, freq_min, freq_max, Nfreq, false, 1.0, B); + // TODO: change here + E = add_dft(direction_component(Ex, d), where->v, freqs, Nfreq, true, 1.0, E); + // TODO: change here + D = add_dft(direction_component(Dx, d), where->v, freqs, Nfreq, false, 1.0, D); + // TODO: change here + H = add_dft(direction_component(Hx, d), where->v, freqs, Nfreq, true, 1.0, H); + // TODO: change here + B = add_dft(direction_component(Bx, d), where->v, freqs, Nfreq, false, 1.0, B); } where = where->next; } delete where_save; - return dft_energy(E, H, D, B, freq_min, freq_max, Nfreq, firstvol); + // TODO: change here + return dft_energy(E, H, D, B, freqs, Nfreq, firstvol); +} + +dft_energy fields::add_dft_energy(const volume_list *where_, double freq_min, double freq_max, + int Nfreq) { + if (Nfreq < 1) abort("Nfreq must be at least 1"); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); + std::vector freqs; + for (int i = 0; i < Nfreq; i++) { + freqs[i] = freq_min + dfreq * i; + } + return add_dft_energy(where_, freqs, Nfreq); } void dft_energy::save_hdf5(h5file *file, const char *dprefix) { @@ -575,6 +683,8 @@ void dft_energy::scale_dfts(complex scale) { } void dft_energy::remove() { + // TODO: add freqs to remove? + //delete[] freqs; while (E) { dft_chunk *nxt = E->next_in_dft; delete E; @@ -613,22 +723,48 @@ direction fields::normal_direction(const volume &where) const { return d; } -dft_flux fields::add_dft_flux(direction d, const volume &where, double freq_min, double freq_max, +// TODO: overload fields::add_dft_flux +dft_flux fields::add_dft_flux(direction d, const volume &where, std::vector freqs, int Nfreq, bool use_symmetry) { if (d == NO_DIRECTION) d = normal_direction(where); volume_list vl(where, direction_component(Sx, d)); - dft_flux flux = add_dft_flux(&vl, freq_min, freq_max, Nfreq, use_symmetry); + // TODO: change here + dft_flux flux = add_dft_flux(&vl, freqs, Nfreq, use_symmetry); flux.normal_direction = d; return flux; } +dft_flux fields::add_dft_flux(direction d, const volume &where, double freq_min, double freq_max, + int Nfreq, bool use_symmetry) { + if (Nfreq < 1) abort("Nfreq must be at least 1"); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); + std::vector freqs]; + for (int i = 0; i < Nfreq; i++) { + freqs[i] = freq_min + dfreq * i; + } + return add_dft_flux(d, where, freqs, Nfreq, use_symmetry); +} + +// TODO: overload fields::add_mode_monitor +dft_flux fields::add_mode_monitor(direction d, const volume &where, std::vector freqs, + int Nfreq) { + // TODO: change here + return add_dft_flux(d, where, freqs, Nfreq, /*use_symmetry=*/false); +} + dft_flux fields::add_mode_monitor(direction d, const volume &where, double freq_min, double freq_max, int Nfreq) { - return add_dft_flux(d, where, freq_min, freq_max, Nfreq, /*use_symmetry=*/false); + if (Nfreq < 1) abort("Nfreq must be at least 1"); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); + std::vector freqs; + for (int i = 0; i < Nfreq; i++) { + freqs[i] = freq_min + dfreq * i; + } + return add_dft_flux(d, where, freqs, Nfreq, /*use_symmetry=*/false); } -dft_flux fields::add_dft_flux_box(const volume &where, double freq_min, double freq_max, - int Nfreq) { +// TODO: overload fields::add_dft_flux_box +dft_flux fields::add_dft_flux_box(const volume &where, std::vector freqs, int Nfreq) { volume_list *faces = 0; LOOP_OVER_DIRECTIONS(where.dim, d) { if (where.in_direction(d) > 0) { @@ -642,28 +778,68 @@ dft_flux fields::add_dft_flux_box(const volume &where, double freq_min, double f } } - dft_flux flux = add_dft_flux(faces, freq_min, freq_max, Nfreq); + // TODO: change here + dft_flux flux = add_dft_flux(faces, freqs, Nfreq); delete faces; return flux; } +dft_flux fields::add_dft_flux_box(const volume &where, double freq_min, double freq_max, + int Nfreq) { + if (Nfreq < 1) abort("Nfreq must be at least 1"); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); + std::vector freqs; + for (int i = 0; i < Nfreq; i++) { + freqs[i] = freq_min + dfreq * i; + } + return add_dft_flux_box(where, freqs, Nfreq); +} + +// TODO: overload fields::add_dft_flux_plane +dft_flux fields::add_dft_flux_plane(const volume &where, std::vector freqs, int Nfreq) { + // TODO: change here + return add_dft_flux(NO_DIRECTION, where, freqs, Nfreq); +} + dft_flux fields::add_dft_flux_plane(const volume &where, double freq_min, double freq_max, int Nfreq) { - return add_dft_flux(NO_DIRECTION, where, freq_min, freq_max, Nfreq); + if (Nfreq < 1) abort("Nfreq must be at least 1"); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); + std::vector freqs; + for (int i = 0; i < Nfreq; i++) { + freqs[i] = freq_min + dfreq * i; + } + return add_dft_flux(NO_DIRECTION, where, freqs, Nfreq); +} + +// TODO: overload dft_fields constructor +dft_fields::dft_fields(dft_chunk *chunks_, std::vector freqs_, int Nfreq_, + const volume &where_) + : where(where_) { + chunks = chunks_; + // TODO: change here + Nfreq = Nfreq_; + freqs.assgin(freqs_.begin(), freqs_.end()); } dft_fields::dft_fields(dft_chunk *chunks_, double freq_min_, double freq_max_, int Nfreq_, const volume &where_) : where(where_) { chunks = chunks_; - freq_min = freq_min_; - dfreq = Nfreq_ <= 1 ? 0.0 : (freq_max_ - freq_min_) / (Nfreq_ - 1); Nfreq = Nfreq_; + double dfreqs = (freq_max_ - freq_min_) / (Nfreq_ - 1); + std::vector freqs_; + for (int i = 0; i < Nfreq_; i++) { + freqs_[i] = freq_min_ + dfreqs * i; + } + freqs.assign(freqs_.begin(), freqs_.end()); } void dft_fields::scale_dfts(cdouble scale) { chunks->scale_dft(scale); } void dft_fields::remove() { + // TODO: add freqs to remove? + // delete[] freqs; while (chunks) { dft_chunk *nxt = chunks->next_in_dft; delete chunks; @@ -671,19 +847,33 @@ void dft_fields::remove() { } } +// TODO: overload fields::add_dft_fields dft_fields fields::add_dft_fields(component *components, int num_components, const volume where, - double freq_min, double freq_max, int Nfreq, bool use_centered_grid) { + std::vector freqs, int Nfreq, bool use_centered_grid) { bool include_dV_and_interp_weights = false; bool sqrt_dV_and_interp_weights = false; // default option from meep.hpp (expose to user?) std::complex extra_weight = 1.0; // default option from meep.hpp (expose to user?) cdouble stored_weight = 1.0; dft_chunk *chunks = 0; for (int nc = 0; nc < num_components; nc++) - chunks = add_dft(components[nc], where, freq_min, freq_max, Nfreq, - include_dV_and_interp_weights, stored_weight, chunks, - sqrt_dV_and_interp_weights,extra_weight,use_centered_grid); + // TODO: change here + chunks = add_dft(components[nc], where, freqs, Nfreq, include_dV_and_interp_weights, + stored_weight, chunks, sqrt_dV_and_interp_weights, extra_weight, + use_centered_grid); + + // TODO: change here + return dft_fields(chunks, freqs, Nfreq, where); +} - return dft_fields(chunks, freq_min, freq_max, Nfreq, where); +dft_fields fields::add_dft_fields(component *components, int num_components, const volume where, + double freq_min, double freq_max, int Nfreq, + bool use_centered_grid) { + double dfreq = (freq_max - freq_min) / (Nfreq - 1); + std::vector freqs; + for (int i = 0; i < Nfreq; i++) { + freqs[i] = freq_min + dfreq * i; + } + return add_dft_fields(components, num_components, where, freqs, Nfreq, use_centered_grid); } /***************************************************************/ diff --git a/src/dft_ldos.cpp b/src/dft_ldos.cpp index add49c9b8..03827b6fd 100644 --- a/src/dft_ldos.cpp +++ b/src/dft_ldos.cpp @@ -22,16 +22,31 @@ using namespace std; namespace meep { -dft_ldos::dft_ldos(double freq_min, double freq_max, int Nfreq) { - if (Nfreq <= 1) { - omega_min = (freq_min + freq_max) * pi; - domega = 0; - Nomega = 1; +// TODO: overload dft_ldos constructor +dft_ldos::dft_ldos(std::vector freqs, int Nfreq) { + if (Nfreq < 1) { abort("Nfreq must be at least 1"); } + else { + Nomega = Nfreq; + std::vector om; + for (int i = 0; i < Nfreq; i++) { om[i] = freqs[i] * 2 * pi; } + omegas.assign(om.begin(), om.end()); } + Fdft = new complex[Nomega]; + Jdft = new complex[Nomega]; + for (int i = 0; i < Nomega; ++i) + Fdft[i] = Jdft[i] = 0.0; + Jsum = 1.0; +} + +// TODO: overload dft_ldos constructor +dft_ldos::dft_ldos(double freq_min, double freq_max, int Nfreq) { + if (Nfreq < 1) { abort("Nfreq must be at least 1"); } else { - omega_min = freq_min * 2 * pi; - domega = (freq_max - freq_min) * 2 * pi / (Nfreq - 1); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); Nomega = Nfreq; + std::vector om; + for (int i = 0; i < Nfreq; i++) { om[i] = (freq_min + dfreq * i) * 2 * pi; } + omegas.assign(om.begin(), om.end()); } Fdft = new complex[Nomega]; Jdft = new complex[Nomega]; @@ -130,8 +145,8 @@ void dft_ldos::update(fields &f) { } } for (int i = 0; i < Nomega; ++i) { - complex Ephase = polar(1.0, (omega_min + i * domega) * f.time()) * scale; - complex Hphase = polar(1.0, (omega_min + i * domega) * (f.time() - f.dt / 2)) * scale; + complex Ephase = polar(1.0, (omegas[i]) * f.time()) * scale; + complex Hphase = polar(1.0, (omegas[i]) * (f.time() - f.dt / 2)) * scale; Fdft[i] += Ephase * EJ + Hphase * HJ; // NOTE: take only 1st time dependence: assumes all sources have same J(t) diff --git a/src/meep.hpp b/src/meep.hpp index 4d95e8289..8bf883b1b 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1009,8 +1009,10 @@ class dft_chunk { void operator-=(const dft_chunk &chunk); // the frequencies to loop_in_chunks - double omega_min, domega; + // TODO: change dft_chunk properties int Nomega; + std::vector omegas; + omegas.reserve(Nomega); component c; // component to DFT (possibly transformed by symmetry) @@ -1081,6 +1083,10 @@ void load_dft_hdf5(dft_chunk *dft_chunks, const char *name, h5file *file, const // dft.cpp (normally created with fields::add_dft_flux) class dft_flux { public: + // TODO: add dft_flux constructor + dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, + std::vector fs, int Nf, const volume &where_, direction normal_direction_, + bool use_symmetry_); dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, double fmin, double fmax, int Nf, const volume &where_, direction normal_direction_, bool use_symmetry_); @@ -1103,8 +1109,10 @@ class dft_flux { void remove(); - double freq_min, dfreq; + // TODO: change dft_flux properties int Nfreq; + std::vector freqs; + freqs.reserve(Nfreq); dft_chunk *E, *H; component cE, cH; volume where; @@ -1115,6 +1123,9 @@ class dft_flux { // dft.cpp (normally created with fields::add_dft_energy) class dft_energy { public: + // TODO: add dft_energy constructor + dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, std::vector fs, + int Nf, const volume &where_); dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double fmin, double fmax, int Nf, const volume &where_); dft_energy(const dft_energy &f); @@ -1140,8 +1151,10 @@ class dft_energy { void remove(); - double freq_min, dfreq; + // TODO: change dft_energy properties int Nfreq; + std::vector freqs; + freqs.reserve(Nfreq); dft_chunk *E, *H, *D, *B; volume where; }; @@ -1149,6 +1162,9 @@ class dft_energy { // stress.cpp (normally created with fields::add_dft_force) class dft_force { public: + // TODO: add dft_force constructor + dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, std::vector fs, + int Nf, const volume &where_); dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double fmin, double fmax, int Nf, const volume &where_); dft_force(const dft_force &f); @@ -1167,8 +1183,10 @@ class dft_force { void remove(); - double freq_min, dfreq; + // TODO: change dft_force properties int Nfreq; + std::vector freqs; + freqs.reserve(Nfreq); dft_chunk *offdiag1, *offdiag2, *diag; volume where; }; @@ -1178,6 +1196,10 @@ class dft_near2far { public: /* fourier tranforms of tangential E and H field components in a medium with the given scalar eps and mu */ + // TODO: add dft_near2far constructor + dft_near2far(dft_chunk *F, std::vector fs, int Nf, double eps, double mu, + const volume &where_, const direction periodic_d_[2], const int periodic_n_[2], + const double periodic_k_[2], const double period_[2]); dft_near2far(dft_chunk *F, double fmin, double fmax, int Nf, double eps, double mu, const volume &where_, const direction periodic_d_[2], const int periodic_n_[2], const double periodic_k_[2], const double period_[2]); @@ -1215,8 +1237,10 @@ class dft_near2far { void remove(); - double freq_min, dfreq; + // TODO: change dft_near2far properties int Nfreq; + std::vector freqs; + freqs.reserve(Nfreq); dft_chunk *F; double eps, mu; volume where; @@ -1225,6 +1249,7 @@ class dft_near2far { double periodic_k[2], period[2]; }; +// dft_ldos.cpp /* Class to compute local-density-of-states spectra: the power spectrum P(omega) of the work done by the sources. Specialized to handle only the case where all sources have the same time dependence, which greatly @@ -1233,10 +1258,13 @@ class dft_near2far { store the Fourier transform per point or per current. */ class dft_ldos { public: + // TODO: add dft_lods constructor + dft_ldos(std::vector freqs, int Nfreq); dft_ldos(double freq_min, double freq_max, int Nfreq); ~dft_ldos() { delete[] Fdft; delete[] Jdft; + // delete[] omegas; } void update(fields &f); // to be called after each timestep @@ -1249,21 +1277,27 @@ class dft_ldos { std::complex *Jdft; // Nomega array of J(t) DFT values double Jsum; // sum of |J| over all points public: - double omega_min, domega; + // TODO: change dft_ldos properties int Nomega; + std::vector omegas; + omegas.reserve(Nomega); }; // dft.cpp (normally created with fields::add_dft_fields) class dft_fields { public: + // TODO: add dft_fields constructor + dft_fields(dft_chunk *chunks, std::vector freqs, int Nfreq, const volume &where); dft_fields(dft_chunk *chunks, double freq_min, double freq_max, int Nfreq, const volume &where); void scale_dfts(std::complex scale); void remove(); - double freq_min, dfreq; + // TODO: change dft_fields properties int Nfreq; + std::vector freqs; + freqs.reserve(Nfreq); dft_chunk *chunks; volume where; }; @@ -1721,29 +1755,58 @@ class fields { double max_abs(derived_component c, const volume &where); // dft.cpp + // TODO: overload add_dft + dft_chunk *add_dft(component c, const volume &where, std::vector freqs, int Nfreq, + bool include_dV_and_interp_weights = true, + std::complex stored_weight = 1.0, dft_chunk *chunk_next = 0, + bool sqrt_dV_and_interp_weights = false, + std::complex extra_weight = 1.0, bool use_centered_grid = true, + int vc = 0); dft_chunk *add_dft(component c, const volume &where, double freq_min, double freq_max, int Nfreq, bool include_dV_and_interp_weights = true, std::complex stored_weight = 1.0, dft_chunk *chunk_next = 0, bool sqrt_dV_and_interp_weights = false, std::complex extra_weight = 1.0, bool use_centered_grid = true, int vc = 0); + // TODO: overload add_dft_pt + dft_chunk *add_dft_pt(component c, const vec &where, std::vector freqs, int Nfreq); dft_chunk *add_dft_pt(component c, const vec &where, double freq_min, double freq_max, int Nfreq); + // TODO: overload add_dft + dft_chunk *add_dft(const volume_list *where, std::vector freqs, int Nfreq, + bool include_dV = true); dft_chunk *add_dft(const volume_list *where, double freq_min, double freq_max, int Nfreq, bool include_dV = true); void update_dfts(); + // TODO: overload add_dft_flux + dft_flux add_dft_flux(const volume_list *where, std::vector freqs, int Nfreq, + bool use_symmetry = true); dft_flux add_dft_flux(const volume_list *where, double freq_min, double freq_max, int Nfreq, bool use_symmetry = true); + // TODO: overload add_dft_flux + dft_flux add_dft_flux(direction d, const volume &where, std::vector freqs, int Nfreq, + bool use_symmetry = true); dft_flux add_dft_flux(direction d, const volume &where, double freq_min, double freq_max, int Nfreq, bool use_symmetry = true); + // TODO: overload add_dft_flux_box + dft_flux add_dft_flux_box(const volume &where, std::vector freqs, int Nfreq); dft_flux add_dft_flux_box(const volume &where, double freq_min, double freq_max, int Nfreq); + // TODO: overload add_dft_flux_plane + dft_flux add_dft_flux_plane(const volume &where, std::vector freqs, int Nfreq); dft_flux add_dft_flux_plane(const volume &where, double freq_min, double freq_max, int Nfreq); // a "mode monitor" is just a dft_flux with symmetry reduction turned off. + // TODO: overload add_mode_monitor + dft_flux add_mode_monitor(direction d, const volume &where, std::vector freqs, int Nfreq); dft_flux add_mode_monitor(direction d, const volume &where, double freq_min, double freq_max, int Nfreq); + + // TODO: overload add_dft_fields + dft_fields add_dft_fields(component *components, int num_components, const volume where, + std::vector freqs, int Nfreq, bool use_centered_grid=true); dft_fields add_dft_fields(component *components, int num_components, const volume where, - double freq_min, double freq_max, int Nfreq, bool use_centered_grid=true); + double freq_min, double freq_max, int Nfreq, + bool use_centered_grid=true); /********************************************************/ /* process_dft_component is an intermediate-level */ @@ -1789,12 +1852,19 @@ class fields { void get_mode_mode_overlap(void *mode1_data, void *mode2_data, dft_flux flux, std::complex overlaps[2]); + // TODO: overload add_dft_energy + dft_energy add_dft_energy(const volume_list *where, std::vector freqs, int Nfreq); dft_energy add_dft_energy(const volume_list *where, double freq_min, double freq_max, int Nfreq); // stress.cpp + // TODO: overload add_dft_force + dft_force add_dft_force(const volume_list *where, std::vector freqs, int Nfreq); dft_force add_dft_force(const volume_list *where, double freq_min, double freq_max, int Nfreq); // near2far.cpp + // TODO: overload add_dft_near2far + dft_near2far add_dft_near2far(const volume_list *where, std::vector freqs, int Nfreq, + int Nperiods = 1); dft_near2far add_dft_near2far(const volume_list *where, double freq_min, double freq_max, int Nfreq, int Nperiods = 1); // monitor.cpp diff --git a/src/mpb.cpp b/src/mpb.cpp index ea1637801..cdc73270d 100644 --- a/src/mpb.cpp +++ b/src/mpb.cpp @@ -751,16 +751,15 @@ void fields::add_eigenmode_source(component c0, const src_time &src, direction d /* similarly, if kpoints is non-null it should point to a */ /* caller-allocated array of size num_bands*num_freqs, which on*/ /* return will be populated by the k-vectors for the modes. */ -/***************************************************************/ +/****************omegas.assign(om.begin(), om.end());***********************************************/ void fields::get_eigenmode_coefficients(dft_flux flux, const volume &eig_vol, int *bands, int num_bands, int parity, double eig_resolution, double eigensolver_tol, std::complex *coeffs, double *vgrp, kpoint_func user_kpoint_func, void *user_kpoint_data, vec *kpoints, vec *kdom_list, double *cscale, direction d) { - double freq_min = flux.freq_min; - double dfreq = flux.dfreq; int num_freqs = flux.Nfreq; + std::vector freqs(flux.freqs); bool match_frequency = true; if (flux.use_symmetry && S.multiplicity() > 1 && parity == 0) @@ -779,7 +778,7 @@ void fields::get_eigenmode_coefficients(dft_flux flux, const volume &eig_vol, in /*- call mpb to compute the eigenmode --------------------------*/ /*--------------------------------------------------------------*/ int band_num = bands[nb]; - double freq = freq_min + nf * dfreq; + double freq = freqs[nf]; double kdom[3]; if (user_kpoint_func) kpoint = user_kpoint_func(freq, band_num, user_kpoint_data); am_now_working_on(MPBTime); diff --git a/src/near2far.cpp b/src/near2far.cpp index 1628d77d7..751a57255 100644 --- a/src/near2far.cpp +++ b/src/near2far.cpp @@ -29,14 +29,34 @@ using namespace std; namespace meep { +// TODO: overload dft_near2far constructor +dft_near2far::dft_near2far(dft_chunk *F_, std::vector fs, int Nf, double eps_, double mu_, + const volume &where_, const direction periodic_d_[2], + const int periodic_n_[2], const double periodic_k_[2], + const double period_[2]) + : Nfreq(Nf), F(F_), eps(eps_), mu(mu_), where(where_) { + // TODO: initialize freqs with fs if Nf >= 1 + if (Nf < 1) abort("Nf must be at least 1"); + freqs.assign(fs.begin(), fs.end()); + for (int i = 0; i < 2; ++i) { + periodic_d[i] = periodic_d_[i]; + periodic_n[i] = periodic_n_[i]; + periodic_k[i] = periodic_k_[i]; + period[i] = period_[i]; + } +} + dft_near2far::dft_near2far(dft_chunk *F_, double fmin, double fmax, int Nf, double eps_, double mu_, const volume &where_, const direction periodic_d_[2], const int periodic_n_[2], const double periodic_k_[2], const double period_[2]) : Nfreq(Nf), F(F_), eps(eps_), mu(mu_), where(where_) { - if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5; - freq_min = fmin; - dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); + // TODO: create fs array with which to initalize freqs + if (Nf < 1) abort("Nf must be at least 1"); + double dfreq = (fmax - fmin) / (Nf - 1); + for (int i = 0; i < Nf; i++) { + freqs[i] = fmin + dfreq * i; + } for (int i = 0; i < 2; ++i) { periodic_d[i] = periodic_d_[i]; periodic_n[i] = periodic_n_[i]; @@ -45,9 +65,11 @@ dft_near2far::dft_near2far(dft_chunk *F_, double fmin, double fmax, int Nf, doub } } +// TODO: change here for updated properties of dft_near2far dft_near2far::dft_near2far(const dft_near2far &f) - : freq_min(f.freq_min), dfreq(f.dfreq), Nfreq(f.Nfreq), F(f.F), eps(f.eps), mu(f.mu), - where(f.where) { + : Nfreq(f.Nfreq), F(f.F), eps(f.eps), mu(f.mu), where(f.where) { + if (Nf < 1) abort("Nf must be at least 1"); + freqs.assign(f.freqs.begin(), f.freqs.end()); for (int i = 0; i < 2; ++i) { periodic_d[i] = f.periodic_d[i]; periodic_n[i] = f.periodic_n[i]; @@ -57,6 +79,8 @@ dft_near2far::dft_near2far(const dft_near2far &f) } void dft_near2far::remove() { + // TODO: add freqs to remove? + // delete[] freqs; while (F) { dft_chunk *nxt = F->next_in_dft; delete F; @@ -329,7 +353,8 @@ void dft_near2far::farfield_lowlevel(std::complex *EH, const vec &x) { #endif for (int i = 0; i < Nfreq; ++i) { std::complex EH6[6]; - double freq = freq_min + i * dfreq; + // TODO: replace freq_min + dfreq * i with freqs[i] + double freq = freqs[i]; size_t idx_dft = 0; LOOP_OVER_IVECS(f->fc->gv, f->is, f->ie, idx) { IVEC_LOOP_LOC(f->fc->gv, x0); @@ -523,7 +548,8 @@ double *dft_near2far::flux(direction df, const volume &where, double resolution) static double approxeq(double a, double b) { return fabs(a - b) < 0.5e-11 * (fabs(a) + fabs(b)); } -dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, double freq_max, +// TODO: overload fields::add_dft_near2far +dft_near2far fields::add_dft_near2far(const volume_list *where, std::vector freqs, int Nfreq, int Nperiods) { dft_chunk *F = 0; /* E and H chunks*/ double eps = 0, mu = 0; @@ -598,14 +624,27 @@ dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, double s = j == 0 ? 1 : -1; /* sign of n x c */ if (is_electric(c)) s = -s; - F = add_dft(c, w->v, freq_min, freq_max, Nfreq, true, s * w->weight, F, false, 1.0, false, - c0); + // TODO: replace freq_min and freq_max with freqs + F = add_dft(c, w->v, freqs, Nfreq, true, s * w->weight, F, false, 1.0, false, c0); } } } - return dft_near2far(F, freq_min, freq_max, Nfreq, eps, mu, everywhere, periodic_d, periodic_n, - periodic_k, period); + // TODO: replace freq_min and freq_max with freqs + return dft_near2far(F, freqs, Nfreq, eps, mu, everywhere, periodic_d, periodic_n, periodic_k, + period); + +} + +dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, double freq_max, + int Nfreq, int Nperiods) { + if (Nfreq < 1) abort("Nfreq must be at least 1"); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); + std::vector freqs; + for (int i = 0; i < Nfreq; i++) { + freqs[i] = freq_min + dfreq * i; + } + return add_dft_near2far(where, freqs, Nfreq, Nperiods); } } // namespace meep diff --git a/src/stress.cpp b/src/stress.cpp index 4dc05eb76..04562c4a8 100644 --- a/src/stress.cpp +++ b/src/stress.cpp @@ -23,13 +23,28 @@ using namespace std; namespace meep { +// TODO: overload dft_force constructor +dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, + std::vector fs, int Nf, const volume &where_) + : where(where_) { + if (Nf < 1) abort("Nf must be at least 1"); + freqs.assign(fs.begin(), fs.end()); + Nfreq = Nf; + offdiag1 = offdiag1_; + offdiag2 = offdiag2_; + diag = diag_; +} dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double fmin, double fmax, int Nf, const volume &where_) : where(where_) { - if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5; - freq_min = fmin; + // TODO: change here + if (Nf < 1) abort("Nf must be at least 1"); Nfreq = Nf; + double dfreq = (fmax - fmin) / (Nf - 1); + for (int i = 0; i < Nf; i++) { + freqs[i] = fmin + dfreq * i; + } dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); offdiag1 = offdiag1_; offdiag2 = offdiag2_; @@ -38,9 +53,9 @@ dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag } dft_force::dft_force(const dft_force &f) : where(f.where) { - freq_min = f.freq_min; + // TODO: change here Nfreq = f.Nfreq; - dfreq = f.dfreq; + freqs.assign(f.freqs.begin(), f.freqs.end()); offdiag1 = f.offdiag1; offdiag2 = f.offdiag2; diag = f.diag; @@ -48,6 +63,8 @@ dft_force::dft_force(const dft_force &f) : where(f.where) { } void dft_force::remove() { + // TODO: add freqs to remove? + // delete[] freqs; while (offdiag1) { dft_chunk *nxt = offdiag1->next_in_dft; delete offdiag1; @@ -132,8 +149,8 @@ void dft_force::scale_dfts(complex scale) { /* note that the components where->c indicate the direction of the force to be computed, so they should be vector components (such as Ex, Ey, ... or Sx, ...) rather than pseudovectors (like Hx, ...). */ -dft_force fields::add_dft_force(const volume_list *where_, double freq_min, double freq_max, - int Nfreq) { +// TODO: overload fields::add_dft_force +dft_force fields::add_dft_force(const volume_list *where_, std::vector freqs, int Nfreq) { dft_chunk *offdiag1 = 0, *offdiag2 = 0, *diag = 0; volume_list *where = S.reduce(where_); @@ -148,28 +165,46 @@ dft_force fields::add_dft_force(const volume_list *where_, double freq_min, doub if (coordinate_mismatch(gv.dim, fd)) abort("coordinate-type mismatch in add_dft_force"); if (fd != nd) { // off-diagaonal stress-tensor terms - offdiag1 = add_dft(direction_component(Ex, fd), where->v, freq_min, freq_max, Nfreq, true, + // TODO: change here + offdiag1 = add_dft(direction_component(Ex, fd), where->v, freqs, Nfreq, true, where->weight, offdiag1); - offdiag2 = add_dft(direction_component(Ex, nd), where->v, freq_min, freq_max, Nfreq, false, + // TODO: change here + offdiag2 = add_dft(direction_component(Ex, nd), where->v, freqs, Nfreq, false, 1.0, offdiag2); - offdiag1 = add_dft(direction_component(Hx, fd), where->v, freq_min, freq_max, Nfreq, true, + // TODO: change here + offdiag1 = add_dft(direction_component(Hx, fd), where->v, freqs, Nfreq, true, where->weight, offdiag1); - offdiag2 = add_dft(direction_component(Hx, nd), where->v, freq_min, freq_max, Nfreq, false, + // TODO: change here + offdiag2 = add_dft(direction_component(Hx, nd), where->v, freqs, Nfreq, false, 1.0, offdiag2); } else // diagonal stress-tensor terms LOOP_OVER_FIELD_DIRECTIONS(gv.dim, d) { complex weight1 = where->weight * (d == fd ? +0.5 : -0.5); - diag = add_dft(direction_component(Ex, d), where->v, freq_min, freq_max, Nfreq, true, 1.0, - diag, true, weight1, false); - diag = add_dft(direction_component(Hx, d), where->v, freq_min, freq_max, Nfreq, true, 1.0, - diag, true, weight1, false); + // TODO: change here + diag = add_dft(direction_component(Ex, d), where->v, freqs, Nfreq, true, 1.0, diag, true, + weight1, false); + // TODO: change here + diag = add_dft(direction_component(Hx, d), where->v, freqs, Nfreq, true, 1.0, diag, true, + weight1, false); } everywhere = everywhere | where->v; } delete where_save; - return dft_force(offdiag1, offdiag2, diag, freq_min, freq_max, Nfreq, everywhere); + // TODO: change here + return dft_force(offdiag1, offdiag2, diag, freqs, Nfreq, everywhere); +} + +dft_force fields::add_dft_force(const volume_list *where_, double freq_min, double freq_max, + int Nfreq) { + if (Nfreq < 1) abort("Nfreq must be at least 1"); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); + std::vector freqs; + for (int i = 0; i < Nfreq; i++) { + freqs[i] = freq_min + dfreq * i; + } + return add_dft_force(where_, freqs, Nfreq); } } // namespace meep