c=3*10^8;
lemdal=1550E-12;
ommiga=lemdal*c;
syms z;
z=90:1:150;
K=exp(-(z-125)/6);
Vin=100*K;
me=9.3E-31;
e=1.6E-19;
chi=45/180*pi;
q0=4.7E3;alpha=1.6E-7;
H=(z-105)./10;
Ne=sqrt(q0/alpha).*exp(1/2.*(1-H-1./cos(chi).*exp(-H)));
epsilong0=8.85E-12;
ommiga_p=sqrt((Ne.'*e^2)./(me*epsilong0));
elta=sqrt(1-ommiga_p.^2./(ommiga.^2+i*ommiga*Vin));
k=imag(elta);
Aref=lemdal*k;
Aabs=sqrt(-1/2.*(1-ommiga_p.^2./(ommiga.^2+Vin.^2))+1/2.*sqrt((1-ommiga_p.^2./(ommiga.^2+Vin.^2)).^2+(Vin./ommiga.*(ommiga_p.^2./(ommiga.^2+Vin.^2))).^2));
fun=@(z)(Aref+Aabs);
Aatt=integral(fun,90,150,'ArrayValued',true);
vpa(Aatt);
[z1,Ne1] = meshgrid(z,Ne);
mesh(z1,Ne1,Aatt)
brighten(0.5)