


{"id":308461,"date":"2023-09-25T16:17:19","date_gmt":"2023-09-25T16:17:19","guid":{"rendered":"\/forum\/forums\/topic\/q-analysis-script-understanding\/"},"modified":"2023-09-25T16:17:19","modified_gmt":"2023-09-25T16:17:19","slug":"q-analysis-script-understanding","status":"closed","type":"topic","link":"https:\/\/innovationspace.ansys.com\/forum\/forums\/topic\/q-analysis-script-understanding\/","title":{"rendered":"Q analysis script understanding"},"content":{"rendered":"<p>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>##############################################<\/p>\n<p># Q analysis<\/p>\n<p># This script calculates the quality factor from<\/p>\n<p># the slope of the decaying envelope.<\/p>\n<p>#<\/p>\n<p># Input properties<\/p>\n<p># make plots: make plots of signal, envelope, slope, spectrum<\/p>\n<p># 1 for yes, 0 for no<\/p>\n<p># number resonances: number of resonances to look for<\/p>\n<p>#<\/p>\n<p># Output properties<\/p>\n<p># f0: vector of resonant frequencies<\/p>\n<p># Q: vector of Q factors of resonances<\/p>\n<p># delta_Q: error in Q factor<\/p>\n<p>#<\/p>\n<p># Tags: resonator high q analysis quality factor<\/p>\n<p>#<\/p>\n<p># Copyright 2012 Lumerical Solutions Inc<\/p>\n<p>##############################################<\/p>\n<p>&nbsp;<\/p>\n<p># simplify input variable names by removing spaces<\/p>\n<p>number_resonances = %number resonances%;<\/p>\n<p>make_plots = %make plots%;<\/p>\n<p>&nbsp;<\/p>\n<p>min_filter_width = 1; # min width of filter in units of resonance FWHM<\/p>\n<p>max_filter_width = 6; # min width of filter in units of resonance FWHM<\/p>\n<p>filter_width_test_points = 20;<\/p>\n<p>zero_pad = 2^16; # fft zero padding<\/p>\n<p># Note fft zero pad should be a power of 2,<\/p>\n<p># and larger gives more resolution in the<\/p>\n<p>#frequency domain.<\/p>\n<p>&nbsp;<\/p>\n<p>for(N=0; (N+1) &lt;= nx*ny*nz; 0) {<\/p>\n<p>N=N+1;} # up to the sum of # of monitors = (nx*ny*nz)<\/p>\n<p>&nbsp;<\/p>\n<p>#################################################<\/p>\n<p># get the monitor data for the first monitor<\/p>\n<p>#################################################<\/p>\n<p>t = getdata(&#8220;t1&#8243;,&#8221;t&#8221;);<\/p>\n<p>field0_t_Ex = pinch(getdata(&#8220;t1&#8243;,&#8221;Ex&#8221;));<\/p>\n<p>field0_t_Ey = pinch(getdata(&#8220;t1&#8243;,&#8221;Ey&#8221;));<\/p>\n<p>field0_t_Ez = pinch(getdata(&#8220;t1&#8243;,&#8221;Ez&#8221;));<\/p>\n<p>&nbsp;<\/p>\n<p>#################################################<\/p>\n<p># do fft to frequency domain for all monitors<\/p>\n<p>#################################################<\/p>\n<p>w = fftw(t,1,zero_pad);<\/p>\n<p>field_w = matrix(length(w),6*N);<\/p>\n<p>for(i=1:N) {<\/p>\n<p>mname = &#8220;t&#8221; + num2str(i);<\/p>\n<p>for(j=1:6) {<\/p>\n<p>if(almostequal(j,1)) { component = &#8220;Ex&#8221;; }<\/p>\n<p>if(almostequal(j,2)) { component = &#8220;Ey&#8221;; }<\/p>\n<p>if(almostequal(j,3)) { component = &#8220;Ez&#8221;; }<\/p>\n<p>if(almostequal(j,4)) { component = &#8220;Hx&#8221;; }<\/p>\n<p>if(almostequal(j,5)) { component = &#8220;Hy&#8221;; }<\/p>\n<p>if(almostequal(j,6)) { component = &#8220;Hz&#8221;; }<\/p>\n<p>if(j &gt; 3.5) { extra_factor = sqrt(mu0\/eps0); }<\/p>\n<p>else { extra_factor = 1; }<\/p>\n<p>if(havedata(mname,component)) {<\/p>\n<p>&nbsp;<\/p>\n<p>field_w(1:length(w),6*(i-1)+j) = 2*extra_factor*( (1:length(w)) &lt;= (length(w)\/2+0.1)) *<\/p>\n<p>fft(pinch(getdata(mname,component)),1,zero_pad);<\/p>\n<p>}<\/p>\n<p>}<\/p>\n<p>}<\/p>\n<p>&nbsp;<\/p>\n<p>&nbsp;<\/p>\n<p>#################################################<\/p>\n<p># find resonant peaks, including all monitors<\/p>\n<p>#################################################<\/p>\n<p>f_spectrum = sum(abs(field_w)^2,2);<\/p>\n<p>&nbsp;<\/p>\n<p>p = findpeaks(f_spectrum,number_resonances);<\/p>\n<p>f0 = w(p)\/(2*pi);<\/p>\n<p>&nbsp;<\/p>\n<p>#################################################<\/p>\n<p># find quality factors<\/p>\n<p>#################################################<\/p>\n<p>&nbsp;<\/p>\n<p># reserve memory for results<\/p>\n<p>peak_spectra = matrix(length(w),number_resonances);<\/p>\n<p>peak_filters2 = matrix(length(w),number_resonances);<\/p>\n<p>&nbsp;<\/p>\n<p># calculate slope of decay using 40-60% of time signal<\/p>\n<p>tp1 = round(0.4*length(t));<\/p>\n<p>tp2 = round(0.6*length(t));<\/p>\n<p>t2 = t(tp1:tp2);<\/p>\n<p>log_field_all = matrix(tp2-tp1+1,number_resonances);<\/p>\n<p>&nbsp;<\/p>\n<p>Q = matrix(number_resonances);<\/p>\n<p>delta_Q = matrix(number_resonances)+1e300;<\/p>\n<p>&nbsp;<\/p>\n<p># loop over each peak<\/p>\n<p>for(i=1:number_resonances) {<\/p>\n<p># find FWHM of peak<\/p>\n<p>peak_val = f_spectrum(p(i));<\/p>\n<p>continue_search = 1;<\/p>\n<p>for(p1=p(i)-1; (p1&gt;1) &amp; continue_search ; 1) {<\/p>\n<p>if(f_spectrum(p1)&lt;=peak_val\/2) {<\/p>\n<p>continue_search = 0;<\/p>\n<p>} else {<\/p>\n<p>p1 = p1-1;<\/p>\n<p>}<\/p>\n<p>}<\/p>\n<p>continue_search = 1;<\/p>\n<p>for(p2=p(i)+1; (p2&lt;length(w))&amp; continue_search; 1) {<\/p>\n<p>if(f_spectrum(p2)&lt;=peak_val\/2) {<\/p>\n<p>continue_search = 0;<\/p>\n<p>} else {<\/p>\n<p>p2 = p2+1;<\/p>\n<p>}<\/p>\n<p>}<\/p>\n<p>if(p1 &lt; 1) { p1 = 1; }<\/p>\n<p>if(p2 &gt; length(w)) { p2 = length(w); }<\/p>\n<p>FWHM = w(p2)-w(p1);<\/p>\n<p>&nbsp;<\/p>\n<p>for(filter_width=linspace(min_filter_width,max_filter_width,filter_width_test_points)) {<\/p>\n<p># calculate the filter for the peak<\/p>\n<p>peak_filter = exp( -0.5*(w-w(p(i)))^2\/(filter_width*FWHM)^2 );<\/p>\n<p>&nbsp;<\/p>\n<p># inverse fft to get data in time domain<\/p>\n<p>field2_t = 0;<\/p>\n<p>for(mcount=1:6*N) {<\/p>\n<p>field2_t = field2_t + abs(invfft(pinch(field_w,2,mcount)*peak_filter))^2;<\/p>\n<p>}<\/p>\n<p>field2_t = field2_t(tp1:tp2);<\/p>\n<p>log_field = log10(abs(field2_t));<\/p>\n<p>&nbsp;<\/p>\n<p># calculate slope and Q from the slope of the decay<\/p>\n<p># estimate error from the slope<\/p>\n<p>slope = (log_field(2:length(t2))-log_field(1:length(t2)-1))\/<\/p>\n<p>(t(2:length(t2))-t(1:length(t2)-1));<\/p>\n<p>slope_mean = sum(slope)\/length(slope);<\/p>\n<p>slope_delta = sqrt( sum((slope-slope_mean)^2)\/length(slope) );<\/p>\n<p>Q_test = -w(p(i))*log10(exp(1))\/(slope_mean);<\/p>\n<p>delta_Q_test = abs(slope_delta\/slope_mean*Q_test);<\/p>\n<p>if(delta_Q_test &lt; delta_Q(i)) {<\/p>\n<p>Q(i) = Q_test;<\/p>\n<p>delta_Q(i) = delta_Q_test;<\/p>\n<p>&nbsp;<\/p>\n<p># collect data for final plot<\/p>\n<p>peak_spectra(1:length(w),i) = f_spectrum * peak_filter^2;<\/p>\n<p>peak_filters2(1:length(w),i) = peak_filter^2;<\/p>\n<p>log_field_all(1:length(t2),i) = log_field;<\/p>\n<p>&nbsp;<\/p>\n<p>}<\/p>\n<p>}<\/p>\n<p># output summary of peak results to script window<\/p>\n<p>?&#8221;Resonance &#8221; + num2str(i) + &#8220;:&#8221;;<\/p>\n<p>?&#8221; frequency = &#8221; + num2str(w(p(i))\/(2*pi)*1e-12) + &#8220;THz, or &#8220;+num2str(2*pi*c\/w(p(i))*1e9)+&#8221; nm&#8221;;<\/p>\n<p>?&#8221; Q = &#8221; + num2str(Q(i)) +&#8221; +\/- &#8221; + num2str(delta_Q(i));<\/p>\n<p>}<\/p>\n<p>&nbsp;<\/p>\n<p>Q_matrix=Q;<\/p>\n<p>&nbsp;<\/p>\n<p>Q = matrixdataset(&#8220;Q&#8221;);<\/p>\n<p>Q.addparameter(&#8220;lambda&#8221;,c\/f0,&#8221;f&#8221;,f0);<\/p>\n<p>Q.addattribute(&#8220;Q&#8221;,Q_matrix);<\/p>\n<p>Q.addattribute(&#8220;dQ&#8221;,delta_Q);<\/p>\n<p>&nbsp;<\/p>\n<p>spectrum = matrixdataset(&#8220;spectrum&#8221;);<\/p>\n<p>spectrum.addparameter(&#8220;lambda&#8221;,c\/(w\/2\/pi),&#8221;f&#8221;,w\/2\/pi);<\/p>\n<p>spectrum.addattribute(&#8220;spectrum&#8221;,f_spectrum);<\/p>\n<p>&nbsp;<\/p>\n<p>#################################################<\/p>\n<p># plot the results<\/p>\n<p>#################################################<\/p>\n<p>if (make_plots) {<\/p>\n<p># plot signal and envelope for first monitor<\/p>\n<p>field_t_Ex = invfft(pinch(field_w,2,1));<\/p>\n<p>field_t_Ex = field_t_Ex(1:length(t));<\/p>\n<p>field_t_Ey = invfft(pinch(field_w,2,2));<\/p>\n<p>field_t_Ey = field_t_Ey(1:length(t));<\/p>\n<p>field_t_Ez = invfft(pinch(field_w,2,3));<\/p>\n<p>field_t_Ez = field_t_Ez(1:length(t));<\/p>\n<p>plot(t*1e15,field0_t_Ex,abs(field_t_Ex),field0_t_Ey,abs(field_t_Ey),field0_t_Ez,abs(field_t_Ez),&#8221;time (fs)&#8221;,&#8221;field envelope&#8221;);<\/p>\n<p>legend(&#8220;field (Ex)&#8221;,&#8221;envelope (Ex)&#8221;,&#8221;field (Ey)&#8221;,&#8221;envelope (Ey)&#8221;,&#8221;field (Ez)&#8221;,&#8221;envelope (Ez)&#8221;);<\/p>\n<p>&nbsp;<\/p>\n<p># plot the slopes of the decaying fields<\/p>\n<p>plot(t2*1e15,log_field_all,&#8221;time (fs)&#8221;,&#8221;log10(|field(t)|)&#8221;,&#8221;Decay for each resonance&#8221;);<\/p>\n<p>&nbsp;<\/p>\n<p># plot spectra<\/p>\n<p>p1 = find(w,0.8*min(w(p)));<\/p>\n<p>p2 = find(w,1.2*max(w(p)));<\/p>\n<p>f = w\/(2*pi);<\/p>\n<p>plot(f(p1:p2)*1e-12,peak_spectra(p1:p2,1:number_resonances)\/max(f_spectrum)<\/p>\n<p>,&#8221;frequency (THz)&#8221;,&#8221;Arbitrary units&#8221;,&#8221;Spectrum of resonances&#8221;);<\/p>\n<p>plot(f(p1:p2)*1e-12,f_spectrum(p1:p2)\/max(f_spectrum),peak_filters2(p1:p2,1:number_resonances)<\/p>\n<p>,&#8221;frequency (THz)&#8221;,&#8221;Arbitrary units&#8221;,&#8221;Spectrum and filters&#8221;);<\/p>\n<p>}<\/p>\n<p>&nbsp;<\/p>\n","protected":false},"template":"","class_list":["post-308461","topic","type-topic","status-closed","hentry","topic-tag-Qfactor-1"],"aioseo_notices":[],"acf":[],"custom_fields":[{"0":{"_bbp_subscription":["291105","3002"],"_bbp_author_ip":["168.143.243.14"]," _bbp_last_reply_id":["0"]," _bbp_likes_count":["0"],"_btv_view_count":["1040"],"_bbp_topic_status":["unanswered"],"_bbp_status":["publish"],"_bbp_topic_id":["308461"],"_bbp_forum_id":["27833"],"_bbp_engagement":["3002","291105"],"_bbp_voice_count":["2"],"_bbp_reply_count":["6"],"_bbp_last_reply_id":["309091"],"_bbp_last_active_id":["309091"],"_bbp_last_active_time":["2023-09-27 12:57:54"]},"test":"md-haziqgraduate-utm-my"}],"_links":{"self":[{"href":"https:\/\/innovationspace.ansys.com\/forum\/wp-json\/wp\/v2\/topics\/308461","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/innovationspace.ansys.com\/forum\/wp-json\/wp\/v2\/topics"}],"about":[{"href":"https:\/\/innovationspace.ansys.com\/forum\/wp-json\/wp\/v2\/types\/topic"}],"version-history":[{"count":0,"href":"https:\/\/innovationspace.ansys.com\/forum\/wp-json\/wp\/v2\/topics\/308461\/revisions"}],"wp:attachment":[{"href":"https:\/\/innovationspace.ansys.com\/forum\/wp-json\/wp\/v2\/media?parent=308461"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}