My question is how to plot the graph for resonant (a.u. vs wavelength), not the (a.u. vs frequency). Which part I need to change. Urgent..Tqs<\/p>\n
##############################################<\/p>\n
# Q analysis<\/p>\n
# This script calculates the quality factor from<\/p>\n
# the slope of the decaying envelope.<\/p>\n
#<\/p>\n
# Input properties<\/p>\n
# make plots: make plots of signal, envelope, slope, spectrum<\/p>\n
# 1 for yes, 0 for no<\/p>\n
# number resonances: number of resonances to look for<\/p>\n
#<\/p>\n
# Output properties<\/p>\n
# f0: vector of resonant frequencies<\/p>\n
# Q: vector of Q factors of resonances<\/p>\n
# delta_Q: error in Q factor<\/p>\n
#<\/p>\n
# Tags: resonator high q analysis quality factor<\/p>\n
#<\/p>\n
# Copyright 2012 Lumerical Solutions Inc<\/p>\n
##############################################<\/p>\n
<\/p>\n
# simplify input variable names by removing spaces<\/p>\n
number_resonances = %number resonances%;<\/p>\n
make_plots = %make plots%;<\/p>\n
<\/p>\n
min_filter_width = 1; # min width of filter in units of resonance FWHM<\/p>\n
max_filter_width = 6; # min width of filter in units of resonance FWHM<\/p>\n
filter_width_test_points = 20;<\/p>\n
zero_pad = 2^16; # fft zero padding<\/p>\n
# Note fft zero pad should be a power of 2,<\/p>\n
# and larger gives more resolution in the<\/p>\n
#frequency domain.<\/p>\n
<\/p>\n
for(N=0; (N+1) <= nx*ny*nz; 0) {<\/p>\n
N=N+1;} # up to the sum of # of monitors = (nx*ny*nz)<\/p>\n
<\/p>\n
#################################################<\/p>\n
# get the monitor data for the first monitor<\/p>\n
#################################################<\/p>\n
t = getdata(“t1″,”t”);<\/p>\n
field0_t_Ex = pinch(getdata(“t1″,”Ex”));<\/p>\n
field0_t_Ey = pinch(getdata(“t1″,”Ey”));<\/p>\n
field0_t_Ez = pinch(getdata(“t1″,”Ez”));<\/p>\n
<\/p>\n
#################################################<\/p>\n
# do fft to frequency domain for all monitors<\/p>\n
#################################################<\/p>\n
w = fftw(t,1,zero_pad);<\/p>\n
field_w = matrix(length(w),6*N);<\/p>\n
for(i=1:N) {<\/p>\n
mname = “t” + num2str(i);<\/p>\n
for(j=1:6) {<\/p>\n
if(almostequal(j,1)) { component = “Ex”; }<\/p>\n
if(almostequal(j,2)) { component = “Ey”; }<\/p>\n
if(almostequal(j,3)) { component = “Ez”; }<\/p>\n
if(almostequal(j,4)) { component = “Hx”; }<\/p>\n
if(almostequal(j,5)) { component = “Hy”; }<\/p>\n
if(almostequal(j,6)) { component = “Hz”; }<\/p>\n
if(j > 3.5) { extra_factor = sqrt(mu0\/eps0); }<\/p>\n
else { extra_factor = 1; }<\/p>\n
if(havedata(mname,component)) {<\/p>\n
<\/p>\n
field_w(1:length(w),6*(i-1)+j) = 2*extra_factor*( (1:length(w)) <= (length(w)\/2+0.1)) *<\/p>\n
fft(pinch(getdata(mname,component)),1,zero_pad);<\/p>\n
}<\/p>\n
}<\/p>\n
}<\/p>\n
<\/p>\n
<\/p>\n
#################################################<\/p>\n
# find resonant peaks, including all monitors<\/p>\n
#################################################<\/p>\n
f_spectrum = sum(abs(field_w)^2,2);<\/p>\n
<\/p>\n
p = findpeaks(f_spectrum,number_resonances);<\/p>\n
f0 = w(p)\/(2*pi);<\/p>\n
<\/p>\n
#################################################<\/p>\n
# find quality factors<\/p>\n
#################################################<\/p>\n
<\/p>\n
# reserve memory for results<\/p>\n
peak_spectra = matrix(length(w),number_resonances);<\/p>\n
peak_filters2 = matrix(length(w),number_resonances);<\/p>\n
<\/p>\n
# calculate slope of decay using 40-60% of time signal<\/p>\n
tp1 = round(0.4*length(t));<\/p>\n
tp2 = round(0.6*length(t));<\/p>\n
t2 = t(tp1:tp2);<\/p>\n
log_field_all = matrix(tp2-tp1+1,number_resonances);<\/p>\n
<\/p>\n
Q = matrix(number_resonances);<\/p>\n
delta_Q = matrix(number_resonances)+1e300;<\/p>\n
<\/p>\n
# loop over each peak<\/p>\n
for(i=1:number_resonances) {<\/p>\n
# find FWHM of peak<\/p>\n
peak_val = f_spectrum(p(i));<\/p>\n
continue_search = 1;<\/p>\n
for(p1=p(i)-1; (p1>1) & continue_search ; 1) {<\/p>\n
if(f_spectrum(p1)<=peak_val\/2) {<\/p>\n
continue_search = 0;<\/p>\n
} else {<\/p>\n
p1 = p1-1;<\/p>\n
}<\/p>\n
}<\/p>\n
continue_search = 1;<\/p>\n
for(p2=p(i)+1; (p2<length(w))& continue_search; 1) {<\/p>\n
if(f_spectrum(p2)<=peak_val\/2) {<\/p>\n
continue_search = 0;<\/p>\n
} else {<\/p>\n
p2 = p2+1;<\/p>\n
}<\/p>\n
}<\/p>\n
if(p1 < 1) { p1 = 1; }<\/p>\n
if(p2 > length(w)) { p2 = length(w); }<\/p>\n
FWHM = w(p2)-w(p1);<\/p>\n
<\/p>\n
for(filter_width=linspace(min_filter_width,max_filter_width,filter_width_test_points)) {<\/p>\n
# calculate the filter for the peak<\/p>\n
peak_filter = exp( -0.5*(w-w(p(i)))^2\/(filter_width*FWHM)^2 );<\/p>\n
<\/p>\n
# inverse fft to get data in time domain<\/p>\n
field2_t = 0;<\/p>\n
for(mcount=1:6*N) {<\/p>\n
field2_t = field2_t + abs(invfft(pinch(field_w,2,mcount)*peak_filter))^2;<\/p>\n
}<\/p>\n
field2_t = field2_t(tp1:tp2);<\/p>\n
log_field = log10(abs(field2_t));<\/p>\n
<\/p>\n
# calculate slope and Q from the slope of the decay<\/p>\n
# estimate error from the slope<\/p>\n
slope = (log_field(2:length(t2))-log_field(1:length(t2)-1))\/<\/p>\n
(t(2:length(t2))-t(1:length(t2)-1));<\/p>\n
slope_mean = sum(slope)\/length(slope);<\/p>\n
slope_delta = sqrt( sum((slope-slope_mean)^2)\/length(slope) );<\/p>\n
Q_test = -w(p(i))*log10(exp(1))\/(slope_mean);<\/p>\n
delta_Q_test = abs(slope_delta\/slope_mean*Q_test);<\/p>\n
if(delta_Q_test < delta_Q(i)) {<\/p>\n
Q(i) = Q_test;<\/p>\n
delta_Q(i) = delta_Q_test;<\/p>\n
<\/p>\n
# collect data for final plot<\/p>\n
peak_spectra(1:length(w),i) = f_spectrum * peak_filter^2;<\/p>\n
peak_filters2(1:length(w),i) = peak_filter^2;<\/p>\n
log_field_all(1:length(t2),i) = log_field;<\/p>\n
<\/p>\n
}<\/p>\n
}<\/p>\n
# output summary of peak results to script window<\/p>\n
?”Resonance ” + num2str(i) + “:”;<\/p>\n
?” frequency = ” + num2str(w(p(i))\/(2*pi)*1e-12) + “THz, or “+num2str(2*pi*c\/w(p(i))*1e9)+” nm”;<\/p>\n
?” Q = ” + num2str(Q(i)) +” +\/- ” + num2str(delta_Q(i));<\/p>\n
}<\/p>\n
<\/p>\n
Q_matrix=Q;<\/p>\n
<\/p>\n
Q = matrixdataset(“Q”);<\/p>\n
Q.addparameter(“lambda”,c\/f0,”f”,f0);<\/p>\n
Q.addattribute(“Q”,Q_matrix);<\/p>\n
Q.addattribute(“dQ”,delta_Q);<\/p>\n
<\/p>\n
spectrum = matrixdataset(“spectrum”);<\/p>\n
spectrum.addparameter(“lambda”,c\/(w\/2\/pi),”f”,w\/2\/pi);<\/p>\n
spectrum.addattribute(“spectrum”,f_spectrum);<\/p>\n
<\/p>\n
#################################################<\/p>\n
# plot the results<\/p>\n
#################################################<\/p>\n
if (make_plots) {<\/p>\n
# plot signal and envelope for first monitor<\/p>\n
field_t_Ex = invfft(pinch(field_w,2,1));<\/p>\n
field_t_Ex = field_t_Ex(1:length(t));<\/p>\n
field_t_Ey = invfft(pinch(field_w,2,2));<\/p>\n
field_t_Ey = field_t_Ey(1:length(t));<\/p>\n
field_t_Ez = invfft(pinch(field_w,2,3));<\/p>\n
field_t_Ez = field_t_Ez(1:length(t));<\/p>\n
plot(t*1e15,field0_t_Ex,abs(field_t_Ex),field0_t_Ey,abs(field_t_Ey),field0_t_Ez,abs(field_t_Ez),”time (fs)”,”field envelope”);<\/p>\n
legend(“field (Ex)”,”envelope (Ex)”,”field (Ey)”,”envelope (Ey)”,”field (Ez)”,”envelope (Ez)”);<\/p>\n
<\/p>\n
# plot the slopes of the decaying fields<\/p>\n
plot(t2*1e15,log_field_all,”time (fs)”,”log10(|field(t)|)”,”Decay for each resonance”);<\/p>\n
<\/p>\n
# plot spectra<\/p>\n
p1 = find(w,0.8*min(w(p)));<\/p>\n
p2 = find(w,1.2*max(w(p)));<\/p>\n
f = w\/(2*pi);<\/p>\n
plot(f(p1:p2)*1e-12,peak_spectra(p1:p2,1:number_resonances)\/max(f_spectrum)<\/p>\n
,”frequency (THz)”,”Arbitrary units”,”Spectrum of resonances”);<\/p>\n
plot(f(p1:p2)*1e-12,f_spectrum(p1:p2)\/max(f_spectrum),peak_filters2(p1:p2,1:number_resonances)<\/p>\n
,”frequency (THz)”,”Arbitrary units”,”Spectrum and filters”);<\/p>\n
}<\/p>\n
<\/p>\n","protected":false},"template":"","class_list":["post-308461","topic","type-topic","status-closed","hentry","topic-tag-Qfactor-1"],"aioseo_notices":[],"aioseo_head":"\n\t\t\n\t\n\t\n\t\n\t\n\t\t\n\t\t\n\t\t\n\t\t\n\t\t\n\t\t\n\t\t\n\t\t\n\t\t\n\t\t\n\t\t\n\t\t