#########################################################
# file: usr_source_from_monitor.lsf
#
# Description:
# 	This script runs simulation 1, extracts field data 
#         from the transmission monitor, then uses that data
#         as a source in simulation 2.
#
# Copyright 2009, Lumerical Solutions, Inc.
#########################################################

###########################################
# Run initial simulation. Result from this simulation will 
# be used as the source for the second simulation.


# load and run initial file
load("usr_source_from_monitor1.fsp");
if (!havedata("transmission","x")) {
  run;
}

# image the electric field over the whole taper, and at the output transmission monitor
visualize(getresult("full_profile","E"));
visualize(getresult("transmission","E"));

# calculate the power transmission
T_sim1 = getresult("transmission","T");
?"Transmission from simulation 1: " +num2str(T_sim1.T);


###########################################
# get the E & H fields to use as a source for the 
# 2nd simulation

# get the field data
E_field=getresult("transmission","E");
H_field=getresult("transmission","H");


# To use E&H data in the import source, the E & H fields
# must be contained in a single dataset
# Alternately, if only using E field data, then the 'E_field' dataset
# can be directly loaded into the Import source.
EM = rectilineardataset("EM fields",E_field.x,E_field.y,E_field.z);
EM.addparameter("lambda",c/E_field.f,"f",E_field.f);   # Note: Import source can only load datasets with one frequency point of data.
EM.addattribute("E",E_field.Ex,E_field.Ey,E_field.Ez);
EM.addattribute("H",H_field.Hx,H_field.Hy,H_field.Hz);


###########################################
# load the 2nd simulation 
load("usr_source_from_monitor2.fsp");


# if not already run, load the monitor data from 
# simulation 1 into the source of simulation 2
if (!havedata("T","x")) {
  select("source2");
  importdataset(EM);

  run;
}


# plot the amount of power absorbed by the particle.  
# Notice that results are normalized by the power transmission 
# from simulation 1.
Pabs = getresult("trans_box","T");
plot(Pabs.lambda*1e9,-Pabs.T*T_sim1.T,"wavelength (nm)","Absorption","Absorption vs wavelength");


# plot the transmission and reflection vs wavelength
T_sim2=getresult("T","T");
R_sim2=getresult("R","T");

plot(c/T_sim2.f*1e9,-T_sim1.T*R_sim2.T,T_sim1.T*T_sim2.T,"wavelength (nm)","R and T","R and T vs wavelength");
legend("R","T");