# file: usr_dipole_power.lsf
# Here we show how a dipole power can be calculated and measured
# in a homogeneous material
# All cases of magnetic and electric dipoles in 2D and 3D are shown
# Copyright (c) 2012 Lumerical Solutions, Inc.

# user-defined settings ################
    two_d = 1; #set to zero to not do 2D
    three_d = 1; #set to zero to not do 3D
    
    index = 1.5; # not free space
    p0 = 1.6e-19*1e-10; # choose an electric dipole amplitude
    m0 = 1e-20; #choose a magnetic dipole amplitude

   lambda1=1e-6;  # specify wavelength range
   lambda2=2e-6;

   dx = 0.067e-6;   # specify mesh size
# end user-defined settings ############

# START 2D section #####################
  if(two_d) {
    for(i=1:4) { #there are 4 cases in 2D
    	#load the template file
	load("usr_dipole_power2D.fsp");
	# make sure we are in cwnorm state
	cwnorm;
	
	#set the refractive index and the dipole amplitudes for each case
	switchtolayout;
	setnamed("FDTD","background index",index);
            setnamed("FDTD","dx",dx);
            setnamed("FDTD","dy",dx);
            setnamed("dipole","wavelength start",lambda1);
            setnamed("dipole","wavelength stop",lambda2);
	if(i==1) { 
	    ?"TM, Electric dipole, broadband source";
	    setnamed("dipole","theta",0);
	    setnamed("dipole","phi",0);
	    setnamed("dipole","dipole type","Electric"); 
                p0_base=p0/getnamed("dipole","base amplitude");
	    setnamed("dipole","amplitude",p0_base); 
	}
	if(i==2) { 
	    ?"TM, Magnetic dipole, broadband source";
	    setnamed("dipole","theta",90);
	    setnamed("dipole","phi",0);
	    setnamed("dipole","dipole type","Magnetic"); 
                m0_base=m0/getnamed("dipole","base amplitude");
	    setnamed("dipole","amplitude",m0_base);
	}
	if(i==3) { 
	    ?"TE, Electric dipole, broadband source";
	    setnamed("dipole","theta",90);
	    setnamed("dipole","phi",0);
	    setnamed("dipole","dipole type","Electric"); 
                p0_base=p0/getnamed("dipole","base amplitude");
	    setnamed("dipole","amplitude",p0_base);
	}
	if(i==4) { 
	    ?"TE, Magnetic dipole, broadband source";
	    setnamed("dipole","theta",0);
	    setnamed("dipole","phi",0);
	    setnamed("dipole","dipole type","Magnetic");
                m0_base=m0/getnamed("dipole","base amplitude");
	    setnamed("dipole","amplitude",m0_base);
	}
	
	# RUN the simulation
	run;
	
	# create frequeny vector
	f = linspace(c/lambda2,c/lambda1,100);

	# calculate power with the built in sourcepower function
	# This function simply evaluates the analytic formula, so 
	# we expect power1 = power2
            power1=sourcepower(f);

	# re-calculate power by manually evaluating analytic formula
	w = 2*pi*f;
	if(i==1) { 
	    power2 = mu0/(4*pi) * p0^2 * w^3 * (pi/2) ; #TM dipole electric
	    caption = "TM Electric dipole of amplitude " + num2str(p0) + " Cm/m";
	}
	if(i==2) { 
	    power2 = index^2 * mu0/(4*pi) * m0^2 * w^3 * (pi/4) / c^2 ; #TM dipole magnetic
	    caption = "TM Magnetic dipole of amplitude " + num2str(m0) + " Am^2/m";
	}
	if(i==3) { 
	    power2 = mu0/(4*pi) * p0^2 * w^3 * (pi/4) ; #TE dipole electric
	    caption = "TE Electric dipole of amplitude " + num2str(p0) + " Cm/m";
	}
	if(i==4) { 
	    power2 = index^2 * mu0/(4*pi) * m0^2 * w^3 * (pi/2) / c^2 ; #TE dipole magnetic
	    caption = "TE Magnetic dipole of amplitude " + num2str(m0) + " Am^2/m";
	}

	# measure the actual power radiated by the dipole with the dipolepower function
	power3 = dipolepower(f);

	# Plot the results, and note the units of W/m because we are in 2D
	plot(c/f*1e6,power1,power2,power3,"wavelength (microns)","Power (Watts/m)",caption);
	legend("sourcepower function","analytic formula","measured power");
	?"maximum difference between measured and analytic is: " + num2str(max( abs(power3-power2)/power2 )*100) + "%";
	# store relative difference
	if(i==1) { diff1 = (power3-power2)/power2; }
	if(i==2) { diff2 = (power3-power2)/power2; }
	if(i==3) { diff3 = (power3-power2)/power2; }
	if(i==4) { diff4 = (power3-power2)/power2; }

    }
    plot(c/f/dx/index,diff1*100,diff2*100,diff3*100,diff4*100,"points per wavelength","difference (%)","2D difference, analytic vs simulation");
    legend("TM electric","TM magnetic","TE electric","TE magnetic");
  }
# end of 2D section ###############

# start 3D section ################ 
  if(three_d) {
    for(i=1:2) { #1 is electric dipole, 2 is magnetic
    	#load the template file
	load("usr_dipole_power3D.fsp");
	# make sure we are in cwnorm state
	cwnorm;

	#set the refractive index and the dipole amplitudes for each case
	switchtolayout;
	setnamed("FDTD","background index",index);
        setnamed("FDTD","dx",dx);
        setnamed("FDTD","dy",dx);
        setnamed("FDTD","dz",dx);
        setnamed("dipole","wavelength start",lambda1);
        setnamed("dipole","wavelength stop",lambda2);
	if(i==1) { 
	    ?msg = "3D, Electric dipole, gaussian source";
	    setnamed("dipole","dipole type","Electric"); 
                p0_base=p0/getnamed("dipole","base amplitude");
	    setnamed("dipole","amplitude",p0_base); 
	} else { 
	    setnamed("dipole","dipole type","Magnetic"); 
	    ?msg = "3D, Magnetic dipole, gaussian source";
                m0_base=m0/getnamed("dipole","base amplitude");
	    setnamed("dipole","amplitude",m0_base); 
	}

	#run the simulation
	run;
	
	# create frequeny vector
	f = linspace(c/lambda2,c/lambda1,100);

	# calculate power with the built in sourcepower function
	# This function simply evaluates the analytic formula, so 
	# we expect power1 = power2
            power1=sourcepower(f);

	# calculate power based on the analytic formula for each case
	w = 2*pi*f;
	if(i==1) {
	    power2 = index*mu0/(4*pi) * p0^2 * w^4 / (3*c) ;
	    caption = "Electric dipole of amplitude " + num2str(p0) + " Cm";
	} else {
	    power2 = index^3*mu0/(4*pi) * m0^2 * w^4  / (3*c^3) ;
	    caption = "Magnetic dipole of amplitude " + num2str(m0) + " Am^2";
	}
	
	# measure the actual power radiated by the dipole with the dipolepower function
	power3 = dipolepower(f);

	# plot the results
	plot(c/f*1e6,power1,power2,power3,"wavelength (microns)","Power (Watts)",caption);
	legend("sourcepower function","analytic formula","measured power");
	?"maximum difference is: " + num2str(max( abs(power3-power2)/power2 )*100) + "%";
	# store the difference
	if(i==1) { diff5 = (power3-power2)/power2; }
	if(i==2) { diff6 = (power3-power2)/power2; }
    }

    plot(c/f/dx/index,diff5*100,diff6*100,"points per wavelength","difference (%)","3D difference, analytic vs simulation");
    legend("electric dipole","magnetic dipole");
  }