clear ; clc;

%讀取電性色散參數檔案

a=textread('Data.Mats','%s','delimiter', 's');

Debyemodel1=str2num(char(a((1))));

Debyemodel2=str2num(char(a((2))));

Debyemodel3=str2num(char(a((3))));

Debyemodel4=str2num(char(a((4))));

Debyemodel5=str2num(char(a((5))));

Drudemodel1=str2num(char(a((6))));

Drudemodel2=str2num(char(a((7))));

Drudemodel3=str2num(char(a((8))));

Drudemodel4=str2num(char(a((9))));

Drudemodel5=str2num(char(a((10))));

Lorentzmodel1=str2num(char(a((11))));

Lorentzmodel2=str2num(char(a((12))));

Lorentzmodel3=str2num(char(a((13))));

Lorentzmodel4=str2num(char(a((14))));

Lorentzmodel5=str2num(char(a((15))));

Pademodel1=str2num(char(a((16))));

Pademodel2=str2num(char(a((17))));

Pademodel3=str2num(char(a((18))));

Pademodel4=str2num(char(a((19))));

Pademodel5=str2num(char(a((20))));

hDebyemodel1=str2num(char(a((21))));

hDebyemodel2=str2num(char(a((22))));

hDebyemodel3=str2num(char(a((23))));

hDebyemodel4=str2num(char(a((24))));

hDebyemodel5=str2num(char(a((25))));

hDrudemodel1=str2num(char(a((26))));

hDrudemodel2=str2num(char(a((27))));

hDrudemodel3=str2num(char(a((28))));

hDrudemodel4=str2num(char(a((29))));

hDrudemodel5=str2num(char(a((30))));

hLorentzmodel1=str2num(char(a((31))));

hLorentzmodel2=str2num(char(a((32))));

hLorentzmodel3=str2num(char(a((33))));

hLorentzmodel4=str2num(char(a((34))));

hLorentzmodel5=str2num(char(a((35))));

Pademodel1=str2num(char(a((36))));

Pademodel2=str2num(char(a((37))));

Pademodel3=str2num(char(a((38))));

Pademodel4=str2num(char(a((39))));

Pademodel5=str2num(char(a((40))));

epsinf1=str2num(char(a((71))));

sigmametal1=str2num(char(a((72))));

epsinf2=str2num(char(a((73))));

sigmametal2=str2num(char(a((74))));

epsinf3=str2num(char(a((75))));

sigmametal3=str2num(char(a((76))));

epsinf4=str2num(char(a((77))));

sigmametal4=str2num(char(a((78))));

epsinf5=str2num(char(a((79))));

sigmametal5=str2num(char(a((80))));

hepsinf1=str2num(char(a((81))));

hsigmametal1=str2num(char(a((82))));

hepsinf2=str2num(char(a((83))));

hsigmametal2=str2num(char(a((84))));

hepsinf3=str2num(char(a((85))));

hsigmametal3=str2num(char(a((86))));

hepsinf4=str2num(char(a((87))));

hsigmametal4=str2num(char(a((88))));

hepsinf5=str2num(char(a((89))));

hsigmametal5=str2num(char(a((90))));

gg=111;

%將各項係數分配入變數

for i=1:Debyemodel1

Detp1(i)=str2num(char(a((gg))));

gamabye1(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Debyemodel2

Detp2(i)=str2num(char(a((gg))));

gamabye2(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Debyemodel3

Detp3(i)=str2num(char(a((gg))));

gamabye3(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Debyemodel4

Detp4(i)=str2num(char(a((gg))));

gamabye4(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Debyemodel5

Detp5(i)=str2num(char(a((gg))));

gamabye5(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Drudemodel1

Wd1(i)=str2num(char(a((gg))));

gamad1(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Drudemodel2

Wd2(i)=str2num(char(a((gg))));

gamad2(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Drudemodel3

Wd3(i)=str2num(char(a((gg))));

gamad3(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Drudemodel4

Wd4(i)=str2num(char(a((gg))));

gamad4(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Drudemodel5

Wd5(i)=str2num(char(a((gg))));

gamad5(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Lorentzmodel1

eps_p1(i)=str2num(char(a((gg))));

Wp1(i)=str2num(char(a((gg+1))));

gamap1(i)=str2num(char(a((gg+2))));

gg=gg+3;

end

for i=1:Lorentzmodel2

eps_p2(i)=str2num(char(a((gg))));

Wp2(i)=str2num(char(a((gg+1))));

gamap2(i)=str2num(char(a((gg+2))));

gg=gg+3;

end

for i=1:Lorentzmodel3

eps_p3(i)=str2num(char(a((gg))));

Wp3(i)=str2num(char(a((gg+1))));

gamap3(i)=str2num(char(a((gg+2))));

gg=gg+3;

end

for i=1:Lorentzmodel4

eps_p4(i)=str2num(char(a((gg))));

Wp4(i)=str2num(char(a((gg+1))));

gamap4(i)=str2num(char(a((gg+2))));

gg=gg+3;

end

for i=1:Lorentzmodel5

eps_p5(i)=str2num(char(a((gg))));

Wp5(i)=str2num(char(a((gg+1))));

gamap5(i)=str2num(char(a((gg+2))));

gg=gg+3;

end

cc=2.99792458e8;

muz=4.0*pi*1.0D-7;

epsz=1.D0/(muz*cc*cc);

lambda=linspace(100,1000,500)*1e-9;

twopic = 1.883651567308853e+09; % twopic=2*pi*c where c is speed of light

omega = twopic.*(lambda.^(-1)); % angular frequency of light (rad/s)

%=================================================

epsilon=epsinf1-sigmametal1./(j*omega*epsz);

%debye

for i=1:Debyemodel1

epsilon=epsilon-Detp1(i)./(1-j*omega*gamabye1(i));

end

%drude

for i=1:Drudemodel1

epsilon=epsilon-Wd1(i)^2./(omega.^2+j*omega*gamad1(i));

end

%lorentz

for i=1:Lorentzmodel1

    epsilon=epsilon-eps_p1(i)^2./(omega.^2+j.*omega*gamap1(i)-Wp1(i)^2);

end

% epsilon_r 對 lambda

figure(1);hold on;

plot(lambda,real(epsilon),'r','linewidth',3);

plot(lambda,imag(epsilon),'b','linewidth',3);

hold off;

title('Au');

xlabel('\lambda (nm)');

ylabel('\epsilon r');

% 反射率頻譜

n=sqrt( (sqrt(real(epsilon).^2+imag(epsilon).^2)+real(epsilon))/2 );

k=sqrt( (sqrt(real(epsilon).^2+imag(epsilon).^2)-real(epsilon))/2 );

r=(n0-n1)./(n0+n1);

figure(1);

plot(lambda,abs(r).^2,'^b','linewidth',1);

title('R');

xlabel('\lambda (nm)');

ylabel('R');