n=32;        %相對折射率
k=0;
fp=fopen('e:\f12.dat','r');
for a0=0:3:90
   a=a0*pi/180;
   k=k+1;
   c=cos(a);
   g=(n*n+c*c-1)^0.5;
   f(k)=((g-c)/(g+c))^2*(1+((c*(g+c)-1)/(c*(g-c)+1))^2)/2;
   fprintf(fp,'%d %f\n',a0,f(k));
end
fclose(fp);
hold on;
plot(0:3:90,f,'r-');