本人要计算BiBO晶体的非共线的相位匹配角，物理模型为一束泵浦光产生满足相位条件会产生两束子光束，其有一定的夹角。求泵浦光在晶体中的相位匹配角度。

本人主方程

clear all

clc

lambdap=0.39 %um

lambdai=0.78

lambdas=0.78

nxp=f(lambdap)

nxi=f(lambdai)

nxs=f(lambdai)

nyp=g(lambdap)

nyi=g(lambdai)

nys=g(lambdas)

nzp=h(lambdap)

nzi=h(lambdai)

nzs=h(lambdas)

alphaphi=0 %3/180*pi

syms alpha beta thetap phip thetas phis thetai phii thetapie

np=k(thetap,phip,lambdap);

ns=l(thetas,phis,lambdas);

ni=l(thetai,phii,lambdai);

eq0=thetap-pi/2-thetapie

eq1=phip-pi/2

eq2=2*np-cos(alpha)*ns-cos(beta)*ni

eq3=sin(alpha)*ns-sin(beta)*ni

eq4=ns*sin(alpha)-sin(alphaphi)

eq5=cos(thetas)+sin(thetapie)*cos(alpha)

eq6=cos(thetai)+sin(thetapie)*cos(beta)

eq7=tan(phii)*sin(beta)-cos(thetapie)*cos(beta)

eq8=tan(phis)*sin(alpha)+cos(alpha)*cos(thetapie)

[alpha,beta,thetap,phip,thetas,phis,thetai,phii,thetapie]=solve(eq0,eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,alpha,beta,thetap,phip,thetas,phis,thetai,phii,thetapie)

调用函数：

function f=f(x)

%BiBO nx

f=sqrt(3.0740+0.0323/(x^2-0.0316)-0.01337*x^2);

end

%nx

function g=g(x)

%BIBO ny

g=sqrt(3.1685+0.0373/(x^2-0.0346)-0.01750*x^2);

end

function h=h(x)

%BIBO nz

h=sqrt(3.6545+0.0511/(x^2-0.0371)-0.0226*x^2);

end

function k=k(theta,phi,lambda)

kx=sin(theta)*cos(phi) ;

ky=sin(theta)*sin(phi);

kz=cos(theta);

nx=f(lambda);

ny=g(lambda);

nz=h(lambda);

a=nx^-2;

b=ny^-2;

c=nz^-2;

B=-(sin(theta)*cos(phi))^2*(b+c)-(sin(theta)*sin(phi))^2*(a+c)-(cos(theta))^2*(a+b);

C=(sin(theta)*cos(phi))^2*b*c+(sin(theta)*sin(phi))^2*a*c+(cos(theta))^2*a*b;

k=sqrt(2)/sqrt(-B-sqrt(B^2-4*C)) %慢光

function l=l(theta,phi,lambda)

kx=sin(theta)*cos(phi) ;

ky=sin(theta)*sin(phi);

kz=cos(theta);

nx=f(lambda);

ny=g(lambda);

nz=h(lambda);

a=nx^-2;

b=ny^-2;

c=nz^-2;

B=-(sin(theta)*cos(phi))^2*(b+c)-(sin(theta)*sin(phi))^2*(a+c)-(cos(theta))^2*(a+b);

C=(sin(theta)*cos(phi))^2*b*c+(sin(theta)*sin(phi))^2*a*c+(cos(theta))^2*a*b;

l=sqrt(2)/sqrt(-B+sqrt(B^2-4*C)) %e2-快光

解除来无解

可能不能用solve