##########################################################################
# Scriptfile: interpChargeInRect.lsf
#
# Description:
# Extract the rectangular grid from the geometry
#
# Inputs: 
# geometryName:    The name of the geometry (rectangle) to interpolate
# mindelta, maxdelta:  max and min spacing (the average will be used if possible
#
# Outputs:
# xgrid
# ygrid
#
# Copyright 2012, Lumerical Solutions, Inc.
##########################################################################


##########################################################################
# Inputs:

mindelta = 5e-9;
maxdelta = 40e-9;
##########################################################################

#
# Get grid. 
#  Inputs: geometryName, mindelta, maxdelta
#  Outputs: xgrid, ygrid
#
?"Creating interpolation grid";
rectGridFromGeom_cdc;

# 
# Get the mesh
#   Outputs: tri, vtx, nvtx, ntri
#
?"Retrieving mesh data";
getTriMesh_cdc;

?"Retrieving electron and hole concentrations";
rhon = getdata("CHARGE","charge","n");
rhop = getdata("CHARGE","charge","p");

?"Extracting results on V[" + contactName + "] = {" + num2str(v(1)) + ", ..., " + num2str(v(nv)) + "}, " + num2str(nv) + " points.";

nbyv = matrix(nvtx,nv);
pbyv = matrix(nvtx,nv);

for (j = 1:nv) {
    nbyv(1:nvtx,j) = pinch(rhon(1:nvtx,1,vindex(j),1)); 
    pbyv(1:nvtx,j) = pinch(rhop(1:nvtx,1,vindex(j),1)); 
}

#
# Interpolate from triangular mesh to rectilinear mesh
#
?"Starting interpolation";
rhoNi = matrix(nxspan,nyspan,nv);
rhoPi = matrix(nxspan,nyspan,nv);

for (j = 1:nv) {
    rhoNi(1:nxspan,1:nyspan,j) = interptri(tri,vtx,nbyv(1:nvtx,j),xgrid,ygrid);
    rhoPi(1:nxspan,1:nyspan,j) = interptri(tri,vtx,pbyv(1:nvtx,j),xgrid,ygrid);
    ?"step " + num2str(j) + "/" + num2str(nv) + " complete";
}

myinf = 1/0;
test_dn = (rhoNi == myinf) or !finite(rhoNi);
test_dp = (rhoPi == myinf) or !finite(rhoPi);

if (sum(test_dn) or sum(test_dp))  {
  ?("WARNING.  'NAN' values found." + endl + 
          "Setting all NAN's to zero." );
  rhoNi(find(test_dn))=0;
  rhoPi(find(test_dp))=0;
}

#
# Save the waveguide data
# Save the to an .ldf file, have at least two points in the third dimension
# since the (n,k) import object reaquires at least two points in each dimension
# 
# the x,y,z transformation axes
x = xgrid;
y = ygrid;
z = [-0.5e-6,0.5e-6];

nx=length(x);
ny=length(y);
nz=length(z);

n_data= matrix(nx,nz,ny,nv);
p_data= matrix(nx,nz,ny,nv);
simOrientation = getnamed("CHARGE simulation region","dimension"); 
x_data = 0;
y_data = 0;
z_data = 0;
if (simOrientation == '2D Z-Normal') {
  for (a = 1:nz) {
    n_data(1:nx,1:ny,a,1:nv) = rhoNi(1:nx,1:ny,1:nv);
    p_data(1:nx,1:ny,a,1:nv) = rhoPi(1:nx,1:ny,1:nv);
  }
  x_data=x;
  y_data=y;
  z_data=z;
} else { if (simOrientation == '2D Y-Normal') {
  for (a = 1:nz) {
    n_data(1:nx,a,1:ny,1:nv) = rhoNi(1:nx,1:ny,1:nv);
    p_data(1:nx,a,1:ny,1:nv) = rhoPi(1:nx,1:ny,1:nv);
  }
  x_data=x;
  y_data=z;
  z_data=y;
} else { if (simOrientation == '2D X-Normal') {
  for (a = 1:nz) {
    n_data(a,1:nx,1:ny,1:nv) = rhoNi(1:nx,1:ny,1:nv);
    p_data(a,1:nx,1:ny,1:nv) = rhoPi(1:nx,1:ny,1:nv);
  }
  x_data=z;
  y_data=x;
  z_data=y;
} } }

filename = "rectCharge_" + geometryName + ".ldf"; 
savedata(filename,x_data,y_data,z_data,n_data,p_data,v);
?"Saved data to file " + filename;