# this code can be used to manually calculate the
# spectrum result from time monitor data

# get E/H time domain fields
fields = getresult("DGTD::monitor","fields");
E = fields.E;
H = fields.H;

# user can specify frequency range of interest
f = linspace(50e12,1000e12,10000);

# calculate czt of time domain data
Exw = czt(pinch(fields.Ex),fields.time,2*pi*f);
Eyw = czt(pinch(fields.Ey),fields.time,2*pi*f);
Ezw = czt(pinch(fields.Ez),fields.time,2*pi*f);
Hxw = czt(pinch(fields.Hx),fields.time,2*pi*f);
Hyw = czt(pinch(fields.Hy),fields.time,2*pi*f);
Hzw = czt(pinch(fields.Hz),fields.time,2*pi*f);

# package results into a dataset
select("DGTD::monitor");
spectrum = rectilineardataset("spectrum",get("x"),get("y"),get("z"));
spectrum.addparameter("lambda",c/f,"f",f);
spectrum.addattribute("E",Exw,Eyw,Ezw);
spectrum.addattribute("H",Hxw,Hyw,Hzw);

# visualize dataset
visualize(spectrum);