closeall;
clear; 
filename = "nanoparticle_dgtd";
N = 5;  # number of points in sweep
W = linspace(0.1,0.6,N);  # vector for width in um
path1 = filedirectory(currentfilename);
path2 = path1 + "/" + filename + "_sweep/";
cd(path2);
Pc = matrix(N,1);
Pt = matrix(N,1);
P = P_flat = P_full = matrix(N,1);
for(i=1:N) {
    load("sweep_" + num2str(i));
    
    E_near = getresult("DGTD::R","fields");
    divergence_angle = getnamed("DGTD::source","divergence angle")/180*pi;  # source div engle in rad
    surf = createsphericalsurface([0,divergence_angle],[0,2*pi],[0,0,1],1,0.05);
    E_far = near2far(E_near,surf);
   
    x = surf.x;
    y = surf.y;
    z = surf.z;
    ele = surf.elements;
   
    E = E_far.E;
    sE = size(E);
    Ne = sE(1);
    E2 = abs(E(1:Ne,1))^2+abs(E(1:Ne,2))^2+abs(E(1:Ne,3))^2;
    E_far.addattribute("E2",E2);
    #visualize(E_far);
    elements = E_far.elements;
    P(i) = quadtri(elements,[x,y,z],E2);
    
    surf_full = createsphericalsurface([0,pi/2],[0,2*pi],[0,0,1],1,0.05);
    E_far_full = near2far(E_near,surf_full);
    E = E_far_full.E;
    sE = size(E);
    Ne = sE(1);
    E2 = abs(E(1:Ne,1))^2+abs(E(1:Ne,2))^2+abs(E(1:Ne,3))^2;
    E_far_full.addattribute("E2",E2);
    x = surf_full.x;
    y = surf_full.y;
    z = surf_full.z;
    elements = E_far_full.elements;
    P_full(i) = quadtri(elements,[x,y,z],E2);
    
    }
    
cd(path1);
load(filename);
plot(W,P/P_full,"width (um)","collected power (normalized)");
legend("");