load('c:\Documents and Settings\Simon\Desktop\ice.txt')

wavlen = ice(:,1);
nr = ice(:,2);
ni = ice(:,3);
m = complex(nr,ni);
thetai = pi/4;

costhetat = (1-((sin(thetai))./m).^2).^(1/2);
Rp = (abs((costhetat-m.*cos(thetai))./(costhetat+m.*cos(thetai)))).^2;
Rs = (abs((cos(thetai)-m.*costhetat)./(cos(thetai)+m.*costhetat))).^2;
R45 = (Rs+Rp)./2;

figure(4);
loglog(wavlen,R45,wavlen,Rs,wavlen,Rp);
xlabel('Wavelength')
ylabel('Reflection')
title('Reflection Coefficient for 45 Degree Incidence')
legend('R45','Rs','Rp')

%% Problem 5
lda = .6328*10^-6;
num = 2*pi/(lda);
nr = 1.333;
ni = 10^-10;
m = complex(nr,ni);
dia = 1:.1:20;

for i=1:1:191
    r = (dia(i)*10^-6)/2;
    ex = Mie(m,r*num);
    qext(i) = ex(1);
    asy(i) = ex(5);
end

asy = asy';
qext = qext';
dia = dia';

figure(8);
plot(dia,asy);
xlabel('Diameter')
ylabel('Asymmetry Parameter')
title('Asymmetry Parameter of 632.8 nm')
figure(9);
plot(dia,qext);
xlabel('Diameter')
ylabel('Qext')
title('Extinction Efficiency of 632.8 nm')
