求高手看看 哪里错了
clear;
clc
global deltk gx
c=2.977e8; %光速c
lamda1=1.064; %泵浦光1的波长,单位为微米
lamda2=1.342; %泵浦光2的波长,单位为微米
lamdaSFG=0.593;%和频光的波长
w1=2*pi*c/(lamda1*1e-6); %泵浦光1的圆频率
w2=2*pi*c/(lamda2*1e-6);
w3=2*pi*c/(lamdaSFG*1e-6);
T=26; %温度单位是摄氏度
E=0:0.5e4:15e5; %外加电压单位v/m
d33=27e-12; %铌酸锂晶体的非线性系数
gama33=30.9e-12; %铌酸锂晶体的电光系数
[Ne1 No1]=index(lamda1,T);
[Ne2 No2]=index(lamda2,T);
[Ne3 No3]=index(lamdaSFG,T);
K1=Ne1*w1/c;
K2=Ne2*w2/c;
K3=Ne3*w3/c;
Deltk=K3-K2-K1;
dx=pi/Deltk; %每个畴的长度
N=round(1e-2/dx);
fx(1)=-1;
for m=2:N
fx(m)=-1*fx(m-1);
end
n1=Ne1-0.5.*fx.*Ne1^3.*gama33.*E;
n2=Ne2-0.5.*fx.*Ne2^3.*gama33.*E;
n3=Ne3-0.5.*fx.*Ne3^3.*gama33.*E;
k1=n1.*w1./c;
k2=n2.*w2./c;
k3=n3.*w3./c;
del=k3-k2-k1;
g=d33.*fx.*sqrt(w1.*w2.*w3./(n1.*n2.*n3))./c;
Inten=5;%泵浦光1的强度
ebsion=8.85e-12;
P=0.5; %I20/I10的值
AA1=sqrt(2*Inten*1e10/(ebsion*c*w1));
AA2=sqrt(2*P*Inten*1e10/(ebsion*c*w2));
for kk=1:100
A1(1)=AA1(kk);
A2(1)=AA2(kk);
A3(1)=0;
for k=1:N
gx=g(k);
deltk=del(k);
options = odeset('RelTol',1e-5,'AbsTol',[1e-5 1e-5 1e-5]);
[X,Y]=ode45(@sfgfu,[(k-1)*dx k*dx],[A1(k) A2(k) A3(k)],options);
A1(k+1) = Y(end,1);
A2(k+1) = Y(end,2);
A3(k+1) = Y(end,3);
%XX(k+1)=X(end);
end
r1(kk)=abs(A1(end))*abs(A1(end))/(A1(1)^2);
r2(kk)=abs(A2(end))*abs(A2(end))/(A2(1)^2);
r3(kk)=w3.*abs(A3(end)).*abs(A3(end))/(w1*A1(1)^2+w2*A2(1)^2);
end
plot(Inten,r1,'-b',Inten,r2,'-.g',Inten,r3,'.r') AA=Inten';
AB=r1';
AC=r2';
AD=r3';
clear;
clc
global deltk gx
c=2.977e8; %光速c
lamda1=1.064; %泵浦光1的波长,单位为微米
lamda2=1.342; %泵浦光2的波长,单位为微米
lamdaSFG=0.593;%和频光的波长
w1=2*pi*c/(lamda1*1e-6); %泵浦光1的圆频率
w2=2*pi*c/(lamda2*1e-6);
w3=2*pi*c/(lamdaSFG*1e-6);
T=26; %温度单位是摄氏度
E=0:0.5e4:15e5; %外加电压单位v/m
d33=27e-12; %铌酸锂晶体的非线性系数
gama33=30.9e-12; %铌酸锂晶体的电光系数
[Ne1 No1]=index(lamda1,T);
[Ne2 No2]=index(lamda2,T);
[Ne3 No3]=index(lamdaSFG,T);
K1=Ne1*w1/c;
K2=Ne2*w2/c;
K3=Ne3*w3/c;
Deltk=K3-K2-K1;
dx=pi/Deltk; %每个畴的长度
N=round(1e-2/dx);
fx(1)=-1;
for m=2:N
fx(m)=-1*fx(m-1);
end
n1=Ne1-0.5.*fx.*Ne1^3.*gama33.*E;
n2=Ne2-0.5.*fx.*Ne2^3.*gama33.*E;
n3=Ne3-0.5.*fx.*Ne3^3.*gama33.*E;
k1=n1.*w1./c;
k2=n2.*w2./c;
k3=n3.*w3./c;
del=k3-k2-k1;
g=d33.*fx.*sqrt(w1.*w2.*w3./(n1.*n2.*n3))./c;
Inten=5;%泵浦光1的强度
ebsion=8.85e-12;
P=0.5; %I20/I10的值
AA1=sqrt(2*Inten*1e10/(ebsion*c*w1));
AA2=sqrt(2*P*Inten*1e10/(ebsion*c*w2));
for kk=1:100
A1(1)=AA1(kk);
A2(1)=AA2(kk);
A3(1)=0;
for k=1:N
gx=g(k);
deltk=del(k);
options = odeset('RelTol',1e-5,'AbsTol',[1e-5 1e-5 1e-5]);
[X,Y]=ode45(@sfgfu,[(k-1)*dx k*dx],[A1(k) A2(k) A3(k)],options);
A1(k+1) = Y(end,1);
A2(k+1) = Y(end,2);
A3(k+1) = Y(end,3);
%XX(k+1)=X(end);
end
r1(kk)=abs(A1(end))*abs(A1(end))/(A1(1)^2);
r2(kk)=abs(A2(end))*abs(A2(end))/(A2(1)^2);
r3(kk)=w3.*abs(A3(end)).*abs(A3(end))/(w1*A1(1)^2+w2*A2(1)^2);
end
plot(Inten,r1,'-b',Inten,r2,'-.g',Inten,r3,'.r') AA=Inten';
AB=r1';
AC=r2';
AD=r3';

