function fig_pfac


A=[
  0.000  1.01300e+03  281.780
  1.436  8.50284e+02  276.334
  2.840  7.13705e+02  269.452
  4.205  5.99064e+02  261.731
  5.528  5.02838e+02  253.277
  6.806  4.22068e+02  244.059
  8.038  3.54272e+02  235.634
  9.232  2.97367e+02  228.744
 10.399  2.49601e+02  225.526
 11.555  2.09508e+02  223.823
 12.704  1.75856e+02  223.142
 13.853  1.47608e+02  223.403
 15.004  1.23898e+02  223.657
 16.156  1.03997e+02  223.919
 17.311  8.72921e+01  224.272
 18.467  7.32705e+01  224.629
 19.626  6.15013e+01  224.960
 20.787  5.16225e+01  225.263
 21.950  4.33305e+01  225.376
 23.114  3.63704e+01  225.768
 24.281  3.05283e+01  226.477
 25.454  2.56246e+01  227.810
 26.636  2.15086e+01  229.563
 27.828  1.80537e+01  231.659
 29.033  1.51538e+01  234.460
 30.254  1.27197e+01  237.714
 31.493  1.06765e+01  241.136
 32.751  8.96160e+00  244.834
 34.029  7.52212e+00  248.874
 35.329  6.31386e+00  252.921
 36.650  5.29968e+00  256.942
 37.993  4.44840e+00  261.213
 39.358  3.73387e+00  265.199
 40.744  3.13410e+00  268.888
 42.149  2.63068e+00  272.404
 43.571  2.20812e+00  275.041
 45.005  1.85344e+00  276.910
 46.447  1.55572e+00  277.982
 47.893  1.30583e+00  277.997
 49.338  1.09608e+00  277.385
 50.779  9.20017e-01  276.267
 52.212  7.72237e-01  274.126
 53.634  6.48195e-01  271.597
 55.042  5.44077e-01  268.771
 56.435  4.56683e-01  265.275
 57.809  3.83327e-01  261.452
 59.163  3.21754e-01  257.543
 60.497  2.70071e-01  253.428
 61.809  2.26691e-01  249.017
 63.099  1.90278e-01  244.674
 64.367  1.59714e-01  240.418
 65.612  1.34059e-01  235.793
 66.834  1.12526e-01  231.249
 68.033  9.44510e-02  226.843
 69.208  7.92796e-02  222.334
 70.361  6.65451e-02  217.753
 71.490  5.58561e-02  213.310
 72.596  4.68841e-02  208.913
 73.679  3.93532e-02  204.192
 74.738  3.30320e-02  199.538
 75.773  2.77261e-02  194.967
 76.783  2.32726e-02  190.010
 77.768  1.95343e-02  185.201
 78.728  1.63966e-02  180.596
 79.664  1.37628e-02  175.969
 80.577  1.15521e-02  171.548
 81.468  9.69655e-03  167.529
 82.339  8.13902e-03  163.888
 83.192  6.83167e-03  160.530
 84.029  5.73431e-03  157.706
 84.853  4.81322e-03  155.373
 85.665  4.04009e-03  153.553
 86.470  3.39114e-03  152.320
 87.270  2.84643e-03  151.531
 88.067  2.38921e-03  151.193
 88.864  2.00544e-03  151.120
 89.660  1.68331e-03  151.328
 90.459  1.41292e-03  151.911
 91.264  1.18597e-03  153.263
 92.076  9.95470e-04  155.041
 92.900  8.35570e-04  157.187
 93.737  7.01354e-04  160.261
 94.592  5.88697e-04  163.853
 95.467  4.94136e-04  167.818
 96.366  4.14764e-04  172.608
 97.293  3.48142e-04  178.311
 98.251  2.92220e-04  184.487
 99.244  2.45282e-04  191.317
100.278  2.05883e-04  199.971
101.360  1.72812e-04  209.434
102.495  1.45054e-04  219.413
103.688  1.21754e-04  231.558
104.947  1.02197e-04  244.077
106.272  8.57814e-05  256.547
107.666  7.20026e-05  269.542
109.125  6.04370e-05  281.011
110.646  5.07291e-05  292.799
112.236  4.25806e-05  306.330
113.916  3.57410e-05  326.767
115.719  3.00000e-05  352.123
];

global Re 
Re = 6.378e6;

ztan = [5e3,10e3,15e3];

slab = str2mat('5 km','10 km','15 km');

figure(1),clf
figure(2),clf

for i = 1:3
  z    = ztan(i):100:100e3;

  p    = exp(interp1(A(:,1),log(A(:,2)),z/1e3));   % [hPa=mbar]
  t    = interp1(A(:,1),A(:,3),z/1e3); 
  n    = 1+77.9593e-6*p./250;

  c    = (Re+ztan(i))*n(1);
  p    = pfac(z,n,c);

  l    = sqrt((Re+z).^2-(Re+ztan(i))^2);
  dl   = diff(l);
  dl   = dl.*p(1:(length(p)-1));
  l    = [0,cumsum(dl)];

  if i==1
   s = '-';
  elseif i == 2
   s = '--';
  else
   s = '-.';
  end

  figure(1)
  plot(p,z/1e3,s)
  hold on

  figure(2)
  plot(l/1e3,p,s)
  hold on
end

figure(1)
xlabel('Prolongation factor [-]','FontSize',12)
ylabel('Altitude [km]','FontSize',12)
legend(slab)
print fig_pfac_z.eps -deps


figure(2)
ylabel('Prolongation factor [-]','FontSize',12)
xlabel('Distance along LOS [km]','FontSize',12)
legend(slab)
print fig_pfac_l.eps -deps




function h = pfac(z,n,c)

global Re;

%if
 i = 2:length(z);
 h(i) = n(i).*sqrt((Re+z(i)).^2-(c/n(1)).^2)./sqrt(((Re+z(i)).*n(i)).^2-c.^2);
 h(1) = sqrt(n(1)/(n(1)+(n(2)-n(1))*(Re+z(1))/(z(2)-z(1))));

%else
% h    = n.*sqrt((Re+z).^2-(c/n(1)).^2)./sqrt(((Re+z).*n).^2-c.^2);
%end