function [val] = compute_radiated_power(p)

a=p(1);
dfrac=p(2);
d=a*dfrac;
np=8;
rot_theta=20;

system(sprintf("nohup mpirun -np %d meep a=%0.2f d=%0.2f rot-theta=%d bragg_outcoupler.ctl > bragg-optimize-a%0.2f-d%0.2f.out",np,a,d,rot_theta,a,d));
system(sprintf("grep farfield: bragg-optimize-a%0.2f-d%0.2f.out |cut -d , -f2- > bragg-optimize-a%0.2f-d%0.2f.dat",a,d,a,d));

# load far-field data from Meep calculation into arrays for each field component
f = dlmread(sprintf("bragg-optimize-a%0.2f-d%0.2f.dat",a,d),',');
Ex = f(:,2); Ey = f(:,3); Ez = f(:,4);
Hx = f(:,5); Hy = f(:,6); Hz = f(:,7);

# conjugate electric field data
Exc = conj(Ex); Eyc = conj(Ey); Ezc = conj(Ez);

# Poynting vector in Cartesian directions
Px = real(Eyc.*Hz - Ezc.*Hy);
Py = real(Ezc.*Hx - Exc.*Hz);
Pz = real(Exc.*Hy - Eyc.*Hx);

# Poynting vector in XZ plane
Pr = sqrt(Px.^2+Pz.^2);

# minimum cone angle (degrees)
ang_min = 70;

# maximum cone angle (degrees)
ang_max = 80;

ang_min = ang_min*pi/180;
ang_max = ang_max*pi/180;

# normalized in-plane flux in range [0,1]
Pnorm = Pr/max(Pr);

# find the array indices for angular range
idx = find((f(:,1) > ang_min) & (f(:,1) < ang_max));

# fraction of in-plane flux concentrated in angular cone
val = sum(Pnorm(idx))/sum(Pnorm);
  
disp(sprintf("power:, %0.2f, %0.2f, %0.6f",a,d,val));
fflush(stdout);
  
return;
