epsilon0 = 8.854*10^-12;
epsilonr = 8.854*10^-12;
c = 3*10^8;
g = 10^-6;
e0 = 1.6*10^-19;
t = 40;
ec = 3;
h = 6.625*10^-34;
vf = 10^6;
ef = 0.05*e0;
ta = 4.5*10^-13;
f = 0.1 : 0.01 : 2.5;
k0 = 2*pi*f;
w = 2*pi*f;
kf = ef/h/vf;
omiga = h*w/ef + i*h*ta^-1/ef;
sigmar = e0^2/h*t*kf/(24*pi)*omiga*unitstep (omiga - 2);
sigmai = e0^2/h*t*kf/(24*pi*pi)*[
   4/omiga - omiga*log[e, abs[4*ec^2/omiga^2 - 4]]];
sigma = sigamr + i*sigmai;
equ = -epsilonr/epsilon0*
  sqrt[k1^2 - k0^2]/sqrt[k1^2 - epsilonr/epsilon0*k0^2]
-[1 + i*sigma*sqrt[k1^2 - k0^2]/(w*epsilon0)]*
  tanh[g*sqrt[k1^2 - epsilonr/epsilon0*k0^2]];
plot[equ, {k1, 2*pi*10^3 + i*2*pi*10^4, 4.2*pi*10^4 + 2*10^4}]
其中f是变化的从0-2.5thz，我不知道怎么表示，然后equ中所有参数都已知，相求出k1的值，跪求大神帮忙，
进而得到折射率k1/k0，然后k1/k0的实部和虚部随频率f的变化的曲线
老师催的我每天很痛苦，我愿散尽所有金币 返回小木虫查看更多