Delta-wave(YY题)

J - Delta-wave
Time Limit:1000MS     Memory Limit:32768KB     64bit IO Format:%I64d & %I64u

Description

A triangle field is numbered with successive integers in the way shown on the picture below. 



The traveller needs to go from the cell with number M to the cell with number N. The traveller is able to enter the cell through cell edges only, he can not travel from cell to cell through vertices. The number of edges the traveller passes makes the length of the traveller's route. 

Write the program to determine the length of the shortest route connecting cells with numbers N and M. 
 

Input

Input contains two integer numbers M and N in the range from 1 to 1000000000 separated with space(s).
 

Output

Output should contain the length of the shortest route.
 

Sample Input

       
6 12
 

Sample Output

       
3
 


题意:额,很容易懂的吧。假设每一个细胞里面都有一个数字,规律就是上面这样,然后一个细胞只能通过细胞壁到相邻的另一个细胞,问给你两个数字,求他们之间最短的距离是闹哪样。

我感觉是一道开脑洞的YY题。上午刚好看到这题想了一下思路,感觉很麻烦啊。。。然后下午去上课的时候自己睡会,脑海里都是这个图形的模拟,然后刷刷刷起来在本子上写下几行代码,然后自己在测数据,发现几处漏洞,继续duang加特技,然后感觉ok了,只等着回去敲。

首先我是这样做的,求出这两个数分别是在第几行。根据每行最后一位数字是一个平方数可以求出来。然后从上面往下面走有一种情况需要满足,那就是细胞壁相邻,但是duang的一下有些不能直接走下来,要先走到同行相邻的细胞才可以。这里就是a&1!=n&1的情况。然后l,r分别表示第一个数在据这行第一个数和最后一个数的距离。l1,r1分别表示第二个数的。。。。然后从上面一行走到下面一行每次需要2步。然后看l,r会不会在l1,r1这个范围内,是的话,就可以通过(b-a-1)*2直接走到。否则就还需要走到最靠近第二个数的位置再在同行内走l-l1+1步或者r-r1+1步。就是这样了。

#include<limits>
#include<queue>
#include<vector>
#include<list>
#include<map>
#include<set>
#include<deque>
#include<stack>
#include<bitset>
#include<algorithm>
#include<functional>
#include<numeric>
#include<utility>
#include<sstream>
#include<iostream>
#include<iomanip>
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<ctime>

#define LL long long
#define eps 1e-8
#define pi acos(-1)
#define INF 0x7fffffff
#define delta 0.98 //模拟退火递增变量
using namespace std;
int main(){
	int n,m,i;
	while (~scanf("%d%d",&n,&m)){
		int a=0,b=0;
		if (n>m){
			n=n+m;
			m=n-m;
			n=n-m;
		}
		for (i=1;;i++){
			if (a==0 && i*i>=n) a=i;
			if (i*i>m){
				b=i;
				break;
			}
		}
		//cout<<a<<" "<<b<<endl;
		if (a==b){
			cout<<m-n<<endl;
			continue;
		}
		int l,r,l1,r1;
		l=n-(a-1)*(a-1);
		r=a*a-n+1;
		l1=m-(b-1)*(b-1);
		r1=b*b-m+1;
		int ans=1;
		if ((a&1)!=(n&1)){
			l--;
			r--;
			ans++;
		}
		if (l<=l1 && r<=r1){
			ans+=(b-a-1)*2;
			if ((b & 1)==(m & 1)) ans++;
		}	
		else{
			if (l1<l) ans+=(b-a-1)*2+l-l1+1; 
			else ans+=(b-a-1)*2+r-r1+1;
		}
		cout<<ans<<endl;
	}
	return 0;
}



""" 3D Forward Simulation with User-Defined Waveforms ================================================= Here we use the module *simpeg.electromagnetics.time_domain* to predict the TDEM response for a trapezoidal waveform. We consider an airborne survey which uses a horizontal coplanar geometry. For this tutorial, we focus on the following: - How to define the transmitters and receivers - How to define more complicated transmitter waveforms - How to define the time-stepping - How to define the survey - How to solve TDEM problems on an OcTree mesh - How to include topography - The units of the conductivity model and resulting data Please note that we have used a coarse mesh and larger time-stepping to shorten the time of the simulation. Proper discretization in space and time is required to simulate the fields at each time channel with sufficient accuracy. """ ######################################################################### # Import Modules # -------------- # from discretize import TreeMesh from discretize.utils import mkvc, refine_tree_xyz, active_from_xyz from simpeg.utils import plot2Ddata from simpeg import maps import simpeg.electromagnetics.time_domain as tdem from tqdm import tqdm from simpeg.electromagnetics.time_domain import sources, receivers import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt import os save_file = False # sphinx_gallery_thumbnail_number = 3 ############################################################### # Defining Topography # ------------------- # # Here we define surface topography as an (N, 3) numpy array. Topography could # also be loaded from a file. Here we define flat topography, however more # complex topographies can be considered. # xx, yy = np.meshgrid(np.linspace(-3000, 3000, 101), np.linspace(-3000, 3000, 101)) zz = np.zeros(np.shape(xx)) topo_xyz = np.c_[mkvc(xx), mkvc(yy), mkvc(zz)] ############################################################### # Defining the Waveform # --------------------- # # Under *simpeg.electromagnetic.time_domain.sources* # there are a multitude of waveforms that can be defined (VTEM, Ramp-off etc...). # Here, we consider a trapezoidal waveform, which consists of a # linear ramp-on followed by a linear ramp-off. For each waveform, it # is important you are cognizant of the off time!!! # # Define a discrete set of times for which your transmitter is 'on'. Here # the waveform is on from -0.002 s to 0 s. waveform_times = np.linspace(-0.002, 0, 21) # For each waveform type, you must define the necessary set of kwargs. # For the trapezoidal waveform we define the ramp on interval, the # ramp-off interval and the off-time. waveform = tdem.sources.TrapezoidWaveform( ramp_on=np.r_[-0.002, -0.001], ramp_off=np.r_[-0.001, 0.0], off_time=0.0 ) # Uncomment to try a quarter sine wave ramp on, followed by a linear ramp-off. # waveform = tdem.sources.QuarterSineRampOnWaveform( # ramp_on=np.r_[-0.002, -0.001], ramp_off=np.r_[-0.001, 0.], off_time=0. # ) # Uncomment to try a custom waveform (just a linear ramp-off). This requires # defining a function for your waveform. # def wave_function(t): # return - t/(np.max(waveform_times) - np.min(waveform_times)) # # waveform = tdem.sources.RawWaveform(waveform_function=wave_function, off_time=0.) # Evaluate the waveform for each on time. waveform_value = [waveform.eval(t) for t in waveform_times] # Plot the waveform fig = plt.figure(figsize=(10, 4)) ax1 = fig.add_subplot(111) ax1.plot(waveform_times, waveform_value, lw=2) ax1.set_xlabel("Times [s]") ax1.set_ylabel("Waveform value") ax1.set_title("Waveform") ##################################################################### # Create Airborne Survey # ---------------------- # # Here we define the survey used in our simulation. For time domain # simulations, we must define the geometry of the source and its waveform. # the receivers, we define their geometry, the type of field they measure and # the time channels at which they measure the field. For this example, # the survey consists of a uniform grid of airborne measurements. # # Observation times for response (time channels) n_times = 10 time_channels = np.logspace(-5, -2.3, n_times) # Defining transmitter locations n_tx =4 xtx, ytx, ztx = np.meshgrid( np.linspace(-50, 50, n_tx), np.linspace(-50, 50, n_tx), [40] ) source_locations = np.c_[mkvc(xtx), mkvc(ytx), mkvc(ztx)] ntx = np.size(xtx) # Define receiver locations xrx, yrx, zrx = np.meshgrid( np.linspace(-50, 50, n_tx), np.linspace(-50, 50, n_tx), [30] ) receiver_locations = np.c_[mkvc(xrx), mkvc(yrx), mkvc(zrx)] source_list = [] # Create empty list to store sources # Each unique location defines a new transmitter for ii in range(ntx): # Here we define receivers that measure the h-field in A/m dbzdt_receiver = tdem.receivers.PointMagneticFluxTimeDerivative( receiver_locations[ii, :], time_channels, "z" ) receivers_list = [ dbzdt_receiver ] # Make a list containing all receivers even if just one # Must define the transmitter properties and associated receivers source = sources.CircularLoop( receivers_list, location=source_locations[ii], # 线圈中心位置 waveform=waveform, # 使用前面定义的波形 radius=5.0, # 设置线圈半径(单位:m) orientation="z", # 线圈法向方向 current=50.0 # 发射电流(A) ) source_list.append(source) survey = tdem.Survey(source_list) ############################################################### # Create OcTree Mesh # ------------------ # # Here we define the OcTree mesh that is used for this example. # We chose to design a coarser mesh to decrease the run time. # When designing a mesh to solve practical time domain problems: # # - Your smallest cell size should be 10%-20% the size of your smallest diffusion distance # - The thickness of your padding needs to be 2-3 times biggest than your largest diffusion distance # - The diffusion distance is ~1260*np.sqrt(rho*t) # # dh = 2.0 # base cell width dom_width = 400.0 # domain width nbc = 2 ** int(np.round(np.log(dom_width / dh) / np.log(2.0))) # num. base cells # Define the base mesh h = [(dh, nbc)] mesh = TreeMesh([h, h, h], x0="CCC") # Mesh refinement based on topography mesh = refine_tree_xyz( mesh, topo_xyz, octree_levels=[0, 0, 1], method="surface", finalize=False ) # Mesh refinement near transmitters and receivers mesh = refine_tree_xyz( mesh, receiver_locations, octree_levels=[ 2,4], method="radial", finalize=False ) # Refine core mesh region xp, yp, zp = np.meshgrid([-80.0, 80.0], [-80.0, 80.0], [-60.0, 0.0]) xyz = np.c_[mkvc(xp), mkvc(yp), mkvc(zp)] mesh = refine_tree_xyz(mesh, xyz, octree_levels=[0, 2,4], method="box", finalize=False) mesh.finalize() ############################################################### # Create Conductivity Model and Mapping for OcTree Mesh # ----------------------------------------------------- # # Here, we define the electrical properties of the Earth as a conductivity # model. The model consists of a conductive block within a more # resistive background. # # Conductivity in S/m air_conductivity = 1e-8 background_conductivity = 2e-3 block_conductivity = 2e0 # Active cells are cells below the surface. ind_active = active_from_xyz(mesh, topo_xyz) if ind_active.sum() == 0: raise ValueError("No active cells detected below topography — check topo or mesh alignment.") model_map = maps.InjectActiveCells(mesh, ind_active, air_conductivity) # Define the model model = background_conductivity * np.ones(ind_active.sum()) ind_block = ( (mesh.gridCC[ind_active, 0] < 10.0) & (mesh.gridCC[ind_active, 0] > -10.0) & (mesh.gridCC[ind_active, 1] < 10.0) & (mesh.gridCC[ind_active, 1] > -10.0) & (mesh.gridCC[ind_active, 2] > -20.0) & (mesh.gridCC[ind_active, 2] < -15.0) ) model[ind_block] = block_conductivity # Plot log-conductivity model mpl.rcParams.update({"font.size": 12}) fig = plt.figure(figsize=(7, 6)) log_model = np.log10(model) plotting_map = maps.InjectActiveCells(mesh, ind_active, np.nan) ax1 = fig.add_axes([0.13, 0.1, 0.6, 0.85]) mesh.plot_slice( plotting_map * log_model, normal="Y", ax=ax1, ind=int(mesh.h[0].size / 2), grid=True, clim=(np.min(log_model), np.max(log_model)), ) ax1.set_title("Conductivity Model at Y = 0 m") ax2 = fig.add_axes([0.75, 0.1, 0.05, 0.85]) norm = mpl.colors.Normalize(vmin=np.min(log_model), vmax=np.max(log_model)) cbar = mpl.colorbar.ColorbarBase( ax2, norm=norm, orientation="vertical", format="$10^{%.1f}$" ) cbar.set_label("Conductivity [S/m]", rotation=270, labelpad=15, size=12) ###################################################### # Define the Time-Stepping # ------------------------ # # Stuff about time-stepping and some rule of thumb # time_steps = [(5e-6, 10), (1e-5, 10), (5e-5, 10), (1e-4, 20), (5e-4, 10), (1e-3, 5) ] ####################################################################### # Simulation: Time-Domain Response # -------------------------------- # # Here we define the formulation for solving Maxwell's equations. Since we are # measuring the time-derivative of the magnetic flux density and working with # a resistivity model, the EB formulation is the most natural. We must also # remember to define the mapping for the conductivity model. # We defined a waveform 'on-time' is from -0.002 s to 0 s. As a result, we need # to set the start time for the simulation to be at -0.002 s. # simulation = tdem.simulation.Simulation3DMagneticFluxDensity( mesh, survey=survey, sigmaMap=model_map, t0=-0.002 ) # Set the time-stepping for the simulation simulation.time_steps = time_steps ######################################################################## # Predict Data and Plot # --------------------- print("开始进行 3D 正演(dpred)...") # SimPEG 的 dpred 实际上是计算每个源的响应,然后拼接起来 # 我们可以手动循环源列表,并用 tqdm 跟踪进度 dpred_list = [] # 使用 tqdm 包装 simulation.survey.source_list 来显示进度 for src in tqdm(simulation.survey.source_list, desc="计算每个发射源的响应"): # 临时创建一个只包含当前源的 Survey 对象 temp_survey = tdem.Survey([src]) # 将 simulation 对象的 survey 临时设置为 temp_survey original_survey = simulation.survey simulation.survey = temp_survey # 预测当前源的数据 # 注意: SimPEG 型的 dpred 方法需要完整的型,即使只计算一个源 # 如果 SimPEG 版本较新,且支持对单个源进行计算,可能更简洁 # 否则,这种临时替换 survey 的方式是通用的 dpred_list.append(simulation.dpred(model)) # 恢复 simulation 对象的 survey simulation.survey = original_survey # 将所有源的预测数据拼接起来 dpred = np.concatenate(dpred_list) print("3D 正演拟完成。") # Predict data for a given model dpred = simulation.dpred(model) # Data were organized by location, then by time channel dpred_plotting = np.reshape(dpred, (n_tx**2, n_times)) # Plot fig = plt.figure(figsize=(10, 4)) # dB/dt at early time v_max = np.max(np.abs(dpred_plotting[:, 0])) ax11 = fig.add_axes([0.05, 0.05, 0.35, 0.9]) plot2Ddata( receiver_locations[:, 0:2], dpred_plotting[:, 0], ax=ax11, ncontour=30, clim=(-v_max, v_max), contourOpts={"cmap": "bwr"}, ) ax11.set_title("dBz/dt at 0.0001 s") ax12 = fig.add_axes([0.42, 0.05, 0.02, 0.9]) norm1 = mpl.colors.Normalize(vmin=-v_max, vmax=v_max) cbar1 = mpl.colorbar.ColorbarBase( ax12, norm=norm1, orientation="vertical", cmap=mpl.cm.bwr ) cbar1.set_label("$T/s$", rotation=270, labelpad=15, size=12) # dB/dt at later times v_max = np.max(np.abs(dpred_plotting[:, -1])) ax21 = fig.add_axes([0.55, 0.05, 0.35, 0.9]) plot2Ddata( receiver_locations[:, 0:2], dpred_plotting[:, -1], ax=ax21, ncontour=30, clim=(-v_max, v_max), contourOpts={"cmap": "bwr"}, ) ax21.set_title("dBz/dt at 0.001 s") ax22 = fig.add_axes([0.92, 0.05, 0.02, 0.9]) norm2 = mpl.colors.Normalize(vmin=-v_max, vmax=v_max) cbar2 = mpl.colorbar.ColorbarBase( ax22, norm=norm2, orientation="vertical", cmap=mpl.cm.bwr ) cbar2.set_label("$T/s$", rotation=270, labelpad=15, size=12) plt.show() ####################################################### # Optional: Export Data # --------------------- for i in range(n_times): fig = plt.figure(figsize=(6, 4)) v_max = np.max(np.abs(dpred_plotting[:, i])) ax = fig.add_axes([0.1, 0.1, 0.75, 0.8]) plot2Ddata( receiver_locations[:, 0:2], dpred_plotting[:, i], ax=ax, ncontour=30, clim=(-v_max, v_max), contourOpts={"cmap": "bwr"}, ) ax.set_title(f"dBz/dt at t = {time_channels[i]:.6f} s") # 加色条 cax = fig.add_axes([0.88, 0.1, 0.03, 0.8]) norm = mpl.colors.Normalize(vmin=-v_max, vmax=v_max) cbar = mpl.colorbar.ColorbarBase(cax, norm=norm, orientation="vertical", cmap=mpl.cm.bwr) cbar.set_label("$T/s$", rotation=270, labelpad=15, size=12) plt.show() # Write the true model, data and topography # if save_file: dir_path = os.path.dirname(tdem.__file__).split(os.path.sep)[:-3] dir_path.extend(["tutorials", "assets", "tdem"]) dir_path = os.path.sep.join(dir_path) + os.path.sep fname = dir_path + "tdem_topo.txt" np.savetxt(fname, np.c_[topo_xyz], fmt="%.4e") # Write data with 2% noise added fname = dir_path + "tdem_data.obs" dpred = dpred + 0.02 * np.abs(dpred) * np.random.randn(len(dpred)) t_vec = np.kron(np.ones(ntx), time_channels) receiver_locations = np.kron(receiver_locations, np.ones((len(time_channels), 1))) np.savetxt(fname, np.c_[receiver_locations, t_vec, dpred], fmt="%.4e") # Plot true model output_model = plotting_map * model output_model[np.isnan(output_model)] = 1e-8 fname = dir_path + "true_model.txt" np.savetxt(fname, output_model, fmt="%.4e") 此代码中,正演拟time_channels和time_steps如何设置参数
10-21
“ %linked with comsol, this program computes the opto-acoustic coupling %coefficients and phononic energy spectra of an <elliptical> photonic %waveguide %note: the spatial characteristics of the excited phononic modes: for large %cross-sectional dimension (micronscale), the difference from the eigenmode %may be more pronounced (due to the densely-distributed eigmodes) %note: the effect of shape may also be more pronounced at sub-wavelength %scales!!! %% tic; clear all; % computing and geometric parameters %Circle: str1='LWRod_d2h_fsbs_si_Rev'; str2='EWRod_d2h_fsbs_si_Rev'; % str1='LWRod_Circle_Aniso_ang4'; % str2='EWRod_Circle_Aniso_fwd_4'; % str1='LWRod_Ellipse_Aniso_ang16_45'; % str2='EWRod_Ellipse_Aniso_fwd_ang16_45'; % str1='LW_ellipse_circle'; % str2='EW_ellipse_circle'; % wavlen=1.55e-6; % dbNr=3.5; % % scz1=0.6977; % Note: in the case of circle, a, b should be taken a value slightly different % from 1.0 to counteract the influence of degeneration (each time the FEM % is run, different result may be obtained) a=sqrt(1)*1.2; b=sqrt(1)*1.2; radius=1.55e-6/4; % radius=0.19375e-06; % scz2=-2*(2*pi*dbNr/(wavlen)*scz1)/(pi/(radius)); % % scz2=0; dxy1=1.0*1.0e-8; dxy2=0.5*1.0e-8; % dxy1=0.1e-8; % dxy2=0.2e-8; freq_opt=1.9355e14; Qf=1000; %mechanical Q-factor num_eigen_ew=0; Num_ew=100; %total number of elastic eigenmodes % optical wave model initialization: md_lw=mphopen(str1); % md_lw=mphload(str1); % md_lw.param.set('ellipse_a',num2str(a)); % md_lw.param.set('ellipse_b',num2str(b)); % md_lw.param.set('dbRadius',horzcat(num2str(radius), ' [m]')); % md_lw.param.set('scz',num2str(scz1)); % md_lw.study('std1').run; % mphsave(md_lw); %normalization of electric field: % X0=-3e-6:(1e-8):3e-6; % Y0=-3e-6:(1e-8):3e-6; X0=-a*radius:dxy1:a*radius; Y0=-b*radius:dxy1:b*radius; [xx,yy]=meshgrid(X0,Y0); Coords=[xx(:),yy(:)]'; Solnum=1; % solnum of optical eigenmode e1 = mphinterp(md_lw,'ewfd.Ex','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e1=reshape(e1,size(xx)); e2 = mphinterp(md_lw,'ewfd.Ey','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e2=reshape(e2,size(xx)); e3 = mphinterp(md_lw,'ewfd.Ez','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e3=reshape(e3,size(xx)); h1 = mphinterp(md_lw,'ewfd.Hx','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); h1=reshape(h1,size(xx)); h2 = mphinterp(md_lw,'ewfd.Hy','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); h2=reshape(h2,size(xx)); h3 = mphinterp(md_lw,'ewfd.Hz','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); h3=reshape(h3,size(xx)); % 在基础代码的这部分替换原有的Pz计算: % 方法一:坡印廷矢量方法 % P^{(i)} = 2∫∫ d²r ẑ·([e^{(i)}]^* × h^{(i)}) % 计算电场和磁场的叉乘 [e]^* × h cross_x = conj(e2).*h3 - conj(e3).*h2; cross_y = conj(e3).*h1 - conj(e1).*h3; cross_z = conj(e1).*h2 - conj(e2).*h1; % 取z分量 ẑ·([e]^* × h) = cross_z integrand = cross_z; % 计算功率积分 PzI = 2 * trapz(Y0, trapz(X0, integrand, 2)); % 验证功率的虚部(理论上应为实数) if abs(imag(PzI)) > 1e-10 * abs(real(PzI)) warning('坡印廷方法:功率积分有显著虚部 %.2e + %.2ei', real(PzI), imag(PzI)); PzI = real(PzI); % 取实部 end Nrf = sqrt(1/PzI); % normalization factor of electric field % compact spatial domain for interpolation of electric fields X0=-(a*radius):dxy1:(a*radius); Y0=-(b*radius):dxy1:(b*radius); % X0=-(a*radius*0.55):dxy:(a*radius*0.55); % Y0=-(b*radius*0.55):dxy:(b*radius*0.55); [xx,yy]=meshgrid(X0,Y0); Coords=[xx(:),yy(:)]'; e1 = mphinterp(md_lw,'ewfd.Ex','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e1=reshape(e1,size(xx)); e2 = mphinterp(md_lw,'ewfd.Ey','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e2=reshape(e2,size(xx)); e3 = mphinterp(md_lw,'ewfd.Ez','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e3=reshape(e3,size(xx)); d1 = mphinterp(md_lw,'ewfd.Dx','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); d1=reshape(d1,size(xx)); d2 = mphinterp(md_lw,'ewfd.Dy','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); d2=reshape(d2,size(xx)); d3 = mphinterp(md_lw,'ewfd.Dz','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); d3=reshape(d3,size(xx)); e1=Nrf*e1; e2=Nrf*e2; e3=Nrf*e3; d1=Nrf*d1; d2=Nrf*d2; d3=Nrf*d3; % % save E_Data e1 e2 e3 d1 d2 d3 xx yy % e1=conj(e1); % e2=conj(e2); % e3=conj(e3); % d1=conj(d1); % d2=conj(d2); % d3=conj(d3); save E1_data e1 xx yy save E2_data e2 xx yy save E3_data e3 xx yy save D1_data d1 xx yy save D2_data d2 xx yy save D3_data d3 xx yy toc %% tic %load the elastic wave model, param setting and run: md_ew=mphopen(str2); % md_ew=mphload(str2); % md_ew.param.set('ellipse_a',num2str(a)); % md_ew.param.set('ellipse_b',num2str(b)); % md_ew.param.set('dbRadius',horzcat(num2str(radius), ' [m]')); % md_ew.param.set('scz',num2str(scz2)); % md_ew.study('std1').run; % str2='EWRod_Circle_Aniso_ang0'; % md_ew=mphopen(str2); % % dth=2*pi/100; % theta=-pi:dth:(pi); % % for i=1:length(theta) % Coords1(:,i)=radius*[cos(theta(i)) sin(theta(i))]'; % end % % SolNum_ew=10; % % mi_kernel = mphinterp(md_ew,'mi_kernel','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',SolNum_ew,'coord',Coords1); % % egd_int=mpheval(md_ew,'egd_int','Complexfun','on','ComplexOut','on','solnum',SolNum_ew); % % Nrf_ew=sqrt(1/abs(egd_int.d1(1,1)/2)); % % mi_kernel=Nrf_ew*mi_kernel; % % figure; hold on % plot(theta, real(mi_kernel), 'r'); % plot(theta, imag(mi_kernel), 'g'); coupcoef_data=zeros(Num_ew,3); % array to save the opto-acoustic coupling coefficients %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% f_count=0; tic % for SolNum_ew=1001:2:2000 % for SolNum_ew=1:Num_ew % for SolNum_ew=1:1 for SolNum_ew=1:2:Num_ew % for SolNum_ew=5 data=mpheval(md_ew,'freq','Complexfun','on','ComplexOut','on','Solnum',SolNum_ew); freq=abs(real(data.d1(:,1))) if abs(freq)>1e8 f_count=f_count+1 rp_coef=mpheval(md_ew,'rp_coef','Complexfun','on','ComplexOut','on','solnum',SolNum_ew); % rp_coef=0; % es_coef=mpheval(md_ew,'es_coef','Complexfun','on','ComplexOut','on','solnum',SolNum_ew); egd_int=mpheval(md_ew,'egd_int','Complexfun','on','ComplexOut','on','solnum',SolNum_ew); es_coef=mpheval(md_ew,'es_coef','Complexfun','on','ComplexOut','on','solnum',SolNum_ew); % es_coef_2=mpheval(md_ew,'es_coef_2','Complexfun','on','ComplexOut','on','solnum',SolNum_ew); % es_coef_3=mpheval(md_ew,'es_coef_3','Complexfun','on','ComplexOut','on','solnum',SolNum_ew); % es_coef_4=mpheval(md_ew,'es_coef_4','Complexfun','on','ComplexOut','on','solnum',SolNum_ew); Nrf_ew=sqrt(1/abs(egd_int.d1(1,1)/2)); Imb=Nrf_ew*rp_coef.d1(1,1) %pay attention to the proportional factor % Imb=0; Ies=Nrf_ew*es_coef.d1(1,1) % Ies_1=Nrf_ew*es_coef_1.d1(1,1) % Ies_2=Nrf_ew*es_coef_2.d1(1,1) % Ies_3=Nrf_ew*es_coef_3.d1(1,1) % Ies_4=Nrf_ew*es_coef_4.d1(1,1) % Ies=0; % coupcoef_data(f_count,:)=[freq Imb Ies]; coupcoef_data(f_count,:)=[freq Imb Ies]; end end toc %generating the phonon spectral data: nf=0; nQf=320; ndf=800; for i=1:f_count freq0=abs(coupcoef_data(i,1)); eg_freq=0; freq1=freq0-nQf*freq0/Qf; freq2=freq0+nQf*freq0/Qf; df=(freq2-freq1)/ndf; for j=1:(ndf+1) freq=freq1+(j-1)*df; eg_freq=0; for k=1:f_count % eg_freq=eg_freq+abs(coupcoef_data(k,2)+coupcoef_data(k,3))^2*freq^4/abs(freq^2-(coupcoef_data(k,1)*(1+sqrt(-1)/(2*Qf)))^2)^2; %derivation??? % % eg_freq=eg_freq+abs(coupcoef_data(k,2))^2*freq^4/abs(freq^2-(coupcoef_data(k,1)*(1+sqrt(-1)/(2*Qf)))^2)^2; eg_freq=eg_freq+abs(coupcoef_data(k,2)+coupcoef_data(k,3))^2*coupcoef_data(k,1)^2/(coupcoef_data(k,1)^2+4*Qf^2*(freq-coupcoef_data(k,1))^2); % coupcoef_data(k,1)^2/(coupcoef_data(k,1)^2+4*Qf^2*(freq-coupcoef_data(k,1))^2); % freq^4/abs(freq^2-(coupcoef_data(k,1)*(1+sqrt(-1)/(2*Qf)))^2)^2; %derivation??? end nf=nf+1; eg_freq_data(nf,1)=abs(freq); eg_freq_data(nf,2)=eg_freq/4; end end % eg_freq_data(:,2)=2*2*pi*freq_opt*Qf*eg_freq_data(:,2); eg_freq_data(:,2)=2*pi*freq_opt*Qf*eg_freq_data(:,2); [a,b]=sort(eg_freq_data(:,1)); eg_freq_data(:,1)=eg_freq_data(b,1); eg_freq_data(:,2)=eg_freq_data(b,2); figure;hold on filename='d2h_fsbs_si_opt1_Rev.xlsx'; xlswrite(filename,eg_freq_data); plot(eg_freq_data(:,1),log10(eg_freq_data(:,2)),'r') % [a,b]=sort(eg_freq_data(:,1)); % figure;hold on % plot(eg_freq_data(b,1),log10(eg_freq_data(b,2)),'r') % save data toc”和“ %linked with comsol, this program computes the opto-acoustic coupling %coefficients and phononic energy spectra of an <elliptical> photonic %waveguide %note: the spatial characteristics of the excited phononic modes: for large %cross-sectional dimension (micronscale), the difference from the eigenmode %may be more pronounced (due to the densely-distributed eigmodes) %note: the effect of shape may also be more pronounced at sub-wavelength %scales!!! %% tic; clear all; % computing and geometric parameters %Circle: str1='LWRod_d2h_fsbs_si_Rev'; str2='EWRod_d2h_fsbs_si_Rev'; % str1='LWRod_Circle_Aniso_ang4'; % str2='EWRod_Circle_Aniso_fwd_4'; % str1='LWRod_Ellipse_Aniso_ang16_45'; % str2='EWRod_Ellipse_Aniso_fwd_ang16_45'; % str1='LW_ellipse_circle'; % str2='EW_ellipse_circle'; % wavlen=1.55e-6; % dbNr=3.5; % % scz1=0.6977; % Note: in the case of circle, a, b should be taken a value slightly different % from 1.0 to counteract the influence of degeneration (each time the FEM % is run, different result may be obtained) a=sqrt(1)*1.2; b=sqrt(1)*1.2; radius=1.55e-6/4; % radius=0.19375e-06; % scz2=-2*(2*pi*dbNr/(wavlen)*scz1)/(pi/(radius)); % % scz2=0; dxy1=1.0*1.0e-8; dxy2=0.5*1.0e-8; % dxy1=0.1e-8; % dxy2=0.2e-8; freq_opt=1.9355e14; Qf=1000; %mechanical Q-factor num_eigen_ew=0; Num_ew=100; %total number of elastic eigenmodes % optical wave model initialization: md_lw=mphopen(str1); % md_lw=mphload(str1); % md_lw.param.set('ellipse_a',num2str(a)); % md_lw.param.set('ellipse_b',num2str(b)); % md_lw.param.set('dbRadius',horzcat(num2str(radius), ' [m]')); % md_lw.param.set('scz',num2str(scz1)); % md_lw.study('std1').run; % mphsave(md_lw); %normalization of electric field: % X0=-3e-6:(1e-8):3e-6; % Y0=-3e-6:(1e-8):3e-6; X0=-a*radius:dxy1:a*radius; Y0=-b*radius:dxy1:b*radius; [xx,yy]=meshgrid(X0,Y0); Coords=[xx(:),yy(:)]'; Solnum=1; % solnum of optical eigenmode e1 = mphinterp(md_lw,'ewfd.Ex','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e1=reshape(e1,size(xx)); e2 = mphinterp(md_lw,'ewfd.Ey','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e2=reshape(e2,size(xx)); e3 = mphinterp(md_lw,'ewfd.Ez','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e3=reshape(e3,size(xx)); h1 = mphinterp(md_lw,'ewfd.Hx','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); h1=reshape(h1,size(xx)); h2 = mphinterp(md_lw,'ewfd.Hy','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); h2=reshape(h2,size(xx)); h3 = mphinterp(md_lw,'ewfd.Hz','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); h3=reshape(h3,size(xx)); % 在基础代码的这部分替换原有的Pz计算: % 方法二:电场旋度方法 % P^{(i)} = 1/(-iωμ₀) ∫∫ d²r { [e^*]·[ẑ×(∇×e)] + [e^*]·[∇×(ẑ×e)] } % 常数定义 omega = 2*pi*freq_opt; % 角频率 [rad/s] mu0 = 4*pi*1e-7; % 真空磁导率 [H/m] % 计算电场梯度 dx = X0(2) - X0(1); % 正确的间距计算 dy = Y0(2) - Y0(1); [dE1_dx, dE1_dy] = gradient(e1, dx, dy); [dE2_dx, dE2_dy] = gradient(e2, dx, dy); [dE3_dx, dE3_dy] = gradient(e3, dx, dy); % 计算ẑ×(∇×e) % ∇×e = (∂Ez/∂y, -∂Ez/∂x, ∂Ey/∂x - ∂Ex/∂y) z_cross_curl_e_x = -(-dE3_dx); % = ∂Ez/∂x z_cross_curl_e_y = -(-dE3_dy); % = ∂Ez/∂y (注意符号) % 计算∇×(ẑ×e) % ẑ×e = (-e2, e1, 0) % ∇×(ẑ×e) = (0, 0, ∂e1/∂x + ∂e2/∂y) curl_z_cross_e_z = dE1_dx + dE2_dy; % 被积函数计算 term1 = conj(e1).*z_cross_curl_e_x + conj(e2).*z_cross_curl_e_y; term2 = conj(e3).*curl_z_cross_e_z; integrand = term1 + term2; % 积分并计算功率 P_integral = trapz(Y0, trapz(X0, integrand, 1)); P_total = real(1/(-1i*omega*mu0) * P_integral); Nrf = sqrt(1/P_total); % normalization factor of electric field % compact spatial domain for interpolation of electric fields X0=-(a*radius):dxy1:(a*radius); Y0=-(b*radius):dxy1:(b*radius); % X0=-(a*radius*0.55):dxy:(a*radius*0.55); % Y0=-(b*radius*0.55):dxy:(b*radius*0.55); [xx,yy]=meshgrid(X0,Y0); Coords=[xx(:),yy(:)]'; e1 = mphinterp(md_lw,'ewfd.Ex','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e1=reshape(e1,size(xx)); e2 = mphinterp(md_lw,'ewfd.Ey','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e2=reshape(e2,size(xx)); e3 = mphinterp(md_lw,'ewfd.Ez','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e3=reshape(e3,size(xx)); d1 = mphinterp(md_lw,'ewfd.Dx','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); d1=reshape(d1,size(xx)); d2 = mphinterp(md_lw,'ewfd.Dy','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); d2=reshape(d2,size(xx)); d3 = mphinterp(md_lw,'ewfd.Dz','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); d3=reshape(d3,size(xx)); e1=Nrf*e1; e2=Nrf*e2; e3=Nrf*e3; d1=Nrf*d1; d2=Nrf*d2; d3=Nrf*d3; % % save E_Data e1 e2 e3 d1 d2 d3 xx yy % e1=conj(e1); % e2=conj(e2); % e3=conj(e3); % d1=conj(d1); % d2=conj(d2); % d3=conj(d3); save E1_data e1 xx yy save E2_data e2 xx yy save E3_data e3 xx yy save D1_data d1 xx yy save D2_data d2 xx yy save D3_data d3 xx yy toc %% tic %load the elastic wave model, param setting and run: md_ew=mphopen(str2); % md_ew=mphload(str2); % md_ew.param.set('ellipse_a',num2str(a)); % md_ew.param.set('ellipse_b',num2str(b)); % md_ew.param.set('dbRadius',horzcat(num2str(radius), ' [m]')); % md_ew.param.set('scz',num2str(scz2)); % md_ew.study('std1').run; % str2='EWRod_Circle_Aniso_ang0'; % md_ew=mphopen(str2); % % dth=2*pi/100; % theta=-pi:dth:(pi); % % for i=1:length(theta) % Coords1(:,i)=radius*[cos(theta(i)) sin(theta(i))]'; % end % % SolNum_ew=10; % % mi_kernel = mphinterp(md_ew,'mi_kernel','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',SolNum_ew,'coord',Coords1); % % egd_int=mpheval(md_ew,'egd_int','Complexfun','on','ComplexOut','on','solnum',SolNum_ew); % % Nrf_ew=sqrt(1/abs(egd_int.d1(1,1)/2)); % % mi_kernel=Nrf_ew*mi_kernel; % % figure; hold on % plot(theta, real(mi_kernel), 'r'); % plot(theta, imag(mi_kernel), 'g'); coupcoef_data=zeros(Num_ew,3); % array to save the opto-acoustic coupling coefficients %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% f_count=0; tic % for SolNum_ew=1001:2:2000 % for SolNum_ew=1:Num_ew % for SolNum_ew=1:1 for SolNum_ew=1:2:Num_ew % for SolNum_ew=5 data=mpheval(md_ew,'freq','Complexfun','on','ComplexOut','on','Solnum',SolNum_ew); freq=abs(real(data.d1(:,1))) if abs(freq)>1e8 f_count=f_count+1 rp_coef=mpheval(md_ew,'rp_coef','Complexfun','on','ComplexOut','on','solnum',SolNum_ew); % rp_coef=0; % es_coef=mpheval(md_ew,'es_coef','Complexfun','on','ComplexOut','on','solnum',SolNum_ew); egd_int=mpheval(md_ew,'egd_int','Complexfun','on','ComplexOut','on','solnum',SolNum_ew); es_coef=mpheval(md_ew,'es_coef','Complexfun','on','ComplexOut','on','solnum',SolNum_ew); % es_coef_2=mpheval(md_ew,'es_coef_2','Complexfun','on','ComplexOut','on','solnum',SolNum_ew); % es_coef_3=mpheval(md_ew,'es_coef_3','Complexfun','on','ComplexOut','on','solnum',SolNum_ew); % es_coef_4=mpheval(md_ew,'es_coef_4','Complexfun','on','ComplexOut','on','solnum',SolNum_ew); Nrf_ew=sqrt(1/abs(egd_int.d1(1,1)/2)); Imb=Nrf_ew*rp_coef.d1(1,1) %pay attention to the proportional factor % Imb=0; Ies=Nrf_ew*es_coef.d1(1,1) % Ies_1=Nrf_ew*es_coef_1.d1(1,1) % Ies_2=Nrf_ew*es_coef_2.d1(1,1) % Ies_3=Nrf_ew*es_coef_3.d1(1,1) % Ies_4=Nrf_ew*es_coef_4.d1(1,1) % Ies=0; % coupcoef_data(f_count,:)=[freq Imb Ies]; coupcoef_data(f_count,:)=[freq Imb Ies]; end end toc %generating the phonon spectral data: nf=0; nQf=320; ndf=800; for i=1:f_count freq0=abs(coupcoef_data(i,1)); eg_freq=0; freq1=freq0-nQf*freq0/Qf; freq2=freq0+nQf*freq0/Qf; df=(freq2-freq1)/ndf; for j=1:(ndf+1) freq=freq1+(j-1)*df; eg_freq=0; for k=1:f_count % eg_freq=eg_freq+abs(coupcoef_data(k,2)+coupcoef_data(k,3))^2*freq^4/abs(freq^2-(coupcoef_data(k,1)*(1+sqrt(-1)/(2*Qf)))^2)^2; %derivation??? % % eg_freq=eg_freq+abs(coupcoef_data(k,2))^2*freq^4/abs(freq^2-(coupcoef_data(k,1)*(1+sqrt(-1)/(2*Qf)))^2)^2; eg_freq=eg_freq+abs(coupcoef_data(k,2)+coupcoef_data(k,3))^2*coupcoef_data(k,1)^2/(coupcoef_data(k,1)^2+4*Qf^2*(freq-coupcoef_data(k,1))^2); % coupcoef_data(k,1)^2/(coupcoef_data(k,1)^2+4*Qf^2*(freq-coupcoef_data(k,1))^2); % freq^4/abs(freq^2-(coupcoef_data(k,1)*(1+sqrt(-1)/(2*Qf)))^2)^2; %derivation??? end nf=nf+1; eg_freq_data(nf,1)=abs(freq); eg_freq_data(nf,2)=eg_freq/4; end end % eg_freq_data(:,2)=2*2*pi*freq_opt*Qf*eg_freq_data(:,2); eg_freq_data(:,2)=2*pi*freq_opt*Qf*eg_freq_data(:,2); [a,b]=sort(eg_freq_data(:,1)); eg_freq_data(:,1)=eg_freq_data(b,1); eg_freq_data(:,2)=eg_freq_data(b,2); figure;hold on filename='d2h_fsbs_si_opt1_Rev.xlsx'; xlswrite(filename,eg_freq_data); plot(eg_freq_data(:,1),log10(eg_freq_data(:,2)),'r') % [a,b]=sort(eg_freq_data(:,1)); % figure;hold on % plot(eg_freq_data(b,1),log10(eg_freq_data(b,2)),'r') % save data toc ”根据上述建议修改后的代码如上,通过Matlab运行程序后的图片坐标值相同,纵坐标数值相差较大,图线的形状相似,请结合相关的物理知识解释一下,其中存在的原因,详细具体的解释一下。
10-11
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值