##########################################################
# Calculate the transmission through a portion of the 
# monitor.  This is done by creating an integration 
# filter.  The filter is set to 1 in the region to integrate
# and set to 0 outside the integration region.
mon_name="transmission";
x=getdata(mon_name,"x");
y=getdata(mon_name,"y");
f=getdata(mon_name,"f");
Pz=getdata(mon_name,"Pz");
Pz=real(Pz);


# define an integration region
# three examples are given, a rectangle, a circle and a polygon
# simply uncomment the line corresponding to the shape you want

#integration_shape = "rectangle";
integration_shape = "circle";
#integration_shape = "polygon";


foundshape = 0;
if (integration_shape=="rectangle") {
  foundshape = 1;
  # a rectangular region
  x_min=-500e-9;	# integration range in m, set by user
  x_max=500e-9;
  y_min=-500e-9;
  y_max=500e-9;
  X=meshgridx(x,y);
  Y=meshgridy(x,y);

  filter = (X<x_max) & (X>x_min);
  filter = filter & (Y<y_max) & (Y>y_min);
}
if(integration_shape=="circle") {
  foundshape = 1;
  # a circular region
  x_center=0e-9;	# integration range in m, set by user
  y_center=0;
  radius=500e-9;
  X=meshgridx(x,y);
  Y=meshgridy(x,y);

  filter = (radius)>=sqrt((X-x_center)^2+(Y-y_center)^2);
}
if(integration_shape=="polygon") {
  foundshape = 1;
  
  # define the vertices
  V = 1e-9*[   0, -800;
            -400, -400;
            -400,  400;
               0,  800;
             400,  400;
             400, -400  ];

  X=meshgridx(x,y);
  Y=meshgridy(x,y);
  filter = matrix(length(x),length(y));
  n = size(V);
  for(i=1:n(1)) {
    if(i==1) { im1 = n(1); } else {im1 = i-1;}
    if(V(im1,2) != V(i,2) ) {
      slope = (V(im1,1) - V(i,1)) / (V(im1,2) - V(i,2));
    } else {
      slope = 0;
    }
    temp =  ( ((V(i,2) <= Y) & (Y < V(im1,2)))| 
               ((V(im1,2) <= Y) & (Y < V(i,2)))    ) &
             (X < slope * (Y - V(i,2)) + V(i,1));
    filter = filter*(!temp) + temp*(!filter);
  }
}
if(!foundshape) {
  ?"Could not recognize shape '" + integration_shape+"'";
  break;
}

# create a copy of the filter for each frequency
filter2= matrix(length(x),length(y),1,length(f));
for (i=1:length(f)) {
  filter2(1:length(x),1:length(y),1,i)=filter;
}


# integrate Pz over the filter region
T_filtered = 0.5* integrate(Pz*filter2,1:2,x,y);
T_filtered = T_filtered/sourcepower(f);

# Plot filter and transmission
image(x*1e9,y*1e9,filter,"x (nm)","y (nm)", integration_shape+" integration region");
plot(f/1e12,T_filtered*100,"frequency (THz)","Power","Power through "+integration_shape+" filter");