Skip to content

Commit

Permalink
compute vgrp even if !match_frequency
Browse files Browse the repository at this point in the history
  • Loading branch information
stevengj committed Nov 9, 2018
1 parent b72c132 commit 013f21a
Showing 1 changed file with 16 additions and 16 deletions.
32 changes: 16 additions & 16 deletions src/mpb.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -387,23 +387,23 @@ void *fields::get_eigenmode(double omega_src,
band_num, G[0][0]*k[0],G[1][1]*k[1],G[2][2]*k[2],
sqrt(eigvals[band_num-1]), num_iters);

if (match_frequency) {
// copy desired single eigenvector into scratch arrays
evectmatrix_resize(&W[0], 1, 0);
evectmatrix_resize(&W[1], 1, 0);
for (int i = 0; i < H.n; ++i)
W[0].data[i] = H.data[H.p-1 + i * H.p];

// compute the group velocity in the kdir direction
maxwell_ucross_op(W[0], W[1], mdata, kdir); // W[1] = (dTheta/dk) W[0]
mpb_real vscratch;
evectmatrix_XtY_diag_real(W[0], W[1], &vgrp, &vscratch);
vgrp /= sqrt(eigvals[band_num - 1]);

// return to original size
evectmatrix_resize(&W[0], band_num, 0);
evectmatrix_resize(&W[1], band_num, 0);
// copy desired single eigenvector into scratch arrays
evectmatrix_resize(&W[0], 1, 0);
evectmatrix_resize(&W[1], 1, 0);
for (int i = 0; i < H.n; ++i)
W[0].data[i] = H.data[H.p-1 + i * H.p];

// compute the group velocity in the kdir direction
maxwell_ucross_op(W[0], W[1], mdata, kdir); // W[1] = (dTheta/dk) W[0]
mpb_real vscratch;
evectmatrix_XtY_diag_real(W[0], W[1], &vgrp, &vscratch);
vgrp /= sqrt(eigvals[band_num - 1]);

// return to original size
evectmatrix_resize(&W[0], band_num, 0);
evectmatrix_resize(&W[1], band_num, 0);

if (match_frequency) {
// update k via Newton step
double dkmatch = (sqrt(eigvals[band_num - 1]) - omega_src) / vgrp;
kmatch = kmatch - dkmatch;
Expand Down

0 comments on commit 013f21a

Please sign in to comment.