% ------------------------------------------------------ */
% MacCready theory by backward induction on Mcready value*/
% This program implements all calculations for           */
% "MacCready theory with uncertain lift and limited      */
%  altitude"                                             */
% The program is a matlab program. Every item is a matrix */
% and dot (.) operators operate element by element.      */
% 
% John H. Cochrane May 4 1999
%    john.cochrane@gsb.uchicago.edu                         */
% converted to matlab Jan 2007 
% ------------------------------------------------------ */


% -------------------------------------------------------------------------
% matlab requires functions to be in separate files. Cut out the following
% and put each in a separate file with its name, i.e. vofg.m, s.m,  etc.
% -----------------------------------------------------------------------


function vl = vofg(m,gl); % returns speed to attain glide g in lift/sink m */
global a smin vms;
    if gl < 0;
        disp('vofg: being asked for negative glide angle');
    end;
    if (m < smin) & (gl > g(m,0));   % can't achieve this glide angle */
        vl = (-1);      % code for can't get there from here */
    else;
        al = a/2;
        bl = -(1/gl+a*vms);
        cl = smin + a/2*vms^2 - m;
        disc = bl^2 - 4*al*cl;
        if disc < 0 ;
            fprintf('vofg error: negative discriminant.  g=%8.3f l=%8.3f gmax= %8.3f\n', gl, m ,G(m,0));
        else;
            vl = (-bl + (disc)^0.5)/(2*al);
        end;
    end;

% -----------------
    
function sink = s(v);
global a vms smin; 
if min(min(v)') < 0;
    fprintf('s(v) error: asking for speed %8,3f less than zero v\n', v);
end;
sink = ( (v >= vms).*(smin+a/2*(v-vms).^2) +...
       (v <  vms).*smin ); %  below min sink, you can circle

% -----------------
   
function ans = rev(x); 
ans = x(end:-1:1,:); 

% -----------------

function yn = mixer(y,x,mix);
% smooth functions across altitude according to normal distribution -- 
% equivalent to adding random altitude per mile */
  if mix == 0;
        yn = y;
  else;
      yn = zeros(size(y,1),1);
      delx = x(2)-x(1); % assumes same grid spacing for all height! */
      maxind = 3*mix; % look 3 sigma on either side  
      indxs = (-maxind:1:maxind)';
      wts = exp(-1/2*((indxs*delx).^2)/(mix^2));
      i = 1;
       while i <= size(x,1);
          test =  (i+indxs > 0) & (i+indxs < size(x,1)); % limit to bottom and top
          %wl = selif(wts,test);
          wl = wts(test); % use only weights on indices between 0 and max x
          wl = wl/sum(wl);
          %il = selif(indxs,test);
          il = indxs(test);
          yn(i) = sum(wl.*y(i+il));
          i = i+1;
      end;
  end;

% -----------------  

function mc_value = mcinv(lam);
global smin a vms; 
    mc_value = ( (lam > -smin).*(abs(2/a*(lam+smin)+vms^2).^0.5)...
          +(lam <= smin).*0);  % speed = 0 or climb in this assumptions */

% -----------------

function mc_value = mc(v);
global a vms smin;
mc_value = (1/2*a*(v.^2 - vms^2) - smin);

% -----------------

function y = interp(x,xgrid,ygrid);
% interpolation proc */
% assumes x is low to high */

    if min(xgrid(2:size(xgrid,1))-xgrid(1:size(xgrid,1)-1)) < 0 ;
        disp('interp: error. xgrid must be organized from low to high');
        end;

    T = size(xgrid,1);
    if x < xgrid(1);
        y = ygrid(1) - (x-xgrid(1))*(ygrid(2)-ygrid(1))/(xgrid(2)-xgrid(1));
    elseif x > xgrid(end);
        y = ygrid(T) + (x-xgrid(T))*(ygrid(T)-ygrid(T-1))/(xgrid(T)-xgrid(T-1));
    else;

     
        indxv = (1:1:size(xgrid,1))'; % note was rows(hgrid), probably a bug in gauss
        dif = x-xgrid;
        low =  dif >= 0; 
        [trash,indxl] = max(indxv.*low); % only want the index, was maxindc in gauss
        indxh = indxl + 1;
        xl = xgrid(indxl);
        yl = ygrid(indxl);
        xh = xgrid(indxh);
        yh = ygrid(indxh);
        y = (yh-yl)/(xh-xl)*(x-xl)+yl;
    
    
    end;

% -----------------

function gl = g(m,lam);
vl = mcinv(lam-m);
gl = vl./(s(vl)-m);

% -----------------

function gl = g2(m,lam);
% g2 returns inverse glide angle (height lost per mile) */
% this can handle negatives (porpoising) */
% speed is restricted to be at least vms so you can't get infinite up */

global a vms smin;
if (size(m,1) > 1) & (size(lam,1) > 1);
    disp('G2 error: one of m, lam may be a vector but not both');
end;
i = 1;
    while i <= size(m,1);
    ml = m(i);
    j = 1;
    while j <= size(lam,1);
        laml = lam(j);
        vl = mcinv(laml-ml);
        if vl == 0;
            vl = vms;
        end;
        if (i == 1) & (j == 1);
            gl = (s(vl)-ml)/vl;
        else;
            gl = [gl;((s(vl)-ml)/vl)];
        end;
        j = j+1;
    end;
    i = i+1;
end;


% --------------------------
% Actual program begins. There should be no matlab code above this line. 
% ------------------------

clear all
close all

global smin a vms; % polar parameters 
% ------------------------------ */
% switches for how program runs  */
% Before you can see the graphs, you will need to run the program with 
% justgrf = 0, and assumptions set to 1, 2, 3, 4, 5, 6, and 7. This will save
% the results for different cases as data files. The final run, or running again
% with justgrf = 1 will then load all the different cases and make the graphs. 
% -----------------------------  */

justgrf = 0;    % allows you to skip everything and play with graph at end */
                % IF you have alreday run once so the results (lamfin) are in
                %   the workspace or saved to files */
                % current graphs require that you run assumptions 2,5,6
                % before making graphs
thermgf = 0;    % plot thermal distribution? */
polarplt = 0;   % plot polar, etc? */
verbose = 0;    % lots of chat */
                               % preprogrammed assumptionss */
                               
assumptions = 7;               % no prepogrammed assumptions 
                               % 1: 0,4 kts, no porpoising, Illinois;
                               %                               dry discus  */
                               % 2: more realistic Illinois; dry discus  */
                               % 3: 1-26, nimbus, realistic Illinois     */
                               % 4: wet discus, uvalde                   */
                               % 5: Ka6 dry illinois      */
                               % 6: illinois with no 1 knot thermals  
                               % 7: illinois ASH25 
                               
rows = inline('size(x,1)');
cols = inline('size(x,2)'); 
                  
kts_to_kmh = 6000*12*2.54/100/1000;
kts_to_ms = kts_to_kmh*1000/3600; 

% -------------------------- */
% Define glider performance  */
% -------------------------- */

smin = 1.15;   % min sink in kts */
vms  = 42;     % speed for min sink in kts */
vhi = 80;      % high speed polar point */
shi = 3;       % high speed sink in kts */
hcap = 0.96;   % used to scale up or down wining speed assumption */

    % preprogrammed glider assumptionss. */
    % 1: dry discus 2: wet Nimbus 3/asw22 3:1-26  4: wet discus 5: KA6

if (assumptions == 0)| (assumptions == 1) | (assumptions == 2) |(assumptions == 6);
    glider = 1;
end;
if assumptions == 3;
    glider = 3;
end;
if assumptions == 4;
    glider = 4;
end;
if assumptions == 5;
    glider = 5;
end;
if assumptions == 7;
    glider = 6;
end;



if glider == 1;
    smin = 1.15;  % min sink in kts */
    vms  = 42;    % speed for min sink in kts */
    vhi = 80;      % high speed polar point */
    shi = 3;      % high speed sink in kts */
    hcap = 0.96;
end;
if glider == 2;
    smin = 0.95;  % min sink in kts */
    vms  = 50;    % speed for min sink in kts */
    vhi = 100;      % high speed polar point */
    shi = 3;      % high speed sink in kts */
    hcap = 0.82;
end;
if glider == 3;
    % 1-26 values 
    smin = 1.7;    % min sink in kts */
    vms  = 33;    % speed for min sink in kts */
    vhi = 64;      % high speed polar point */
    shi = 4;      % high speed sink in kts */
    hcap = 1.6;    
end;
if glider == 4;
    ldg = (9/6.5)^0.5; % wingloading factor */
    smin = 1.15*ldg;  % min sink in kts */
    vms  = 42*ldg;    % speed for min sink in kts */
    vhi = 80*ldg;      % high speed polar point */
    shi = 3*ldg;      % high speed sink in kts */
    hcap = 1;        % reset V target for this assumptions */
end;

if glider == 5;
    % Ka6CR values 
    smin = 0.7/kts_to_ms;    % min sink in m/s, convert to kts */
    vms  = 65/kts_to_kmh;    % speed for min sink in km/h */
    vhi = 120/kts_to_kmh;      % high speed polar point */
    shi = 1.7/kts_to_ms;      % high speed sink in ms */
    hcap = 1.3;    
end; 

if glider == 6;
    % ASH25 values 
    smin = 0.5/kts_to_ms;    % min sink in m/s, convert to kts */
    vms  = 80/kts_to_kmh;    % speed for min sink in km/h */
    vhi = 150/kts_to_kmh;      % high speed polar point */
    shi = 1/kts_to_ms;      % high speed sink in ms */
    hcap = 0.82;    
end; 


% ------------------------------ */
% Define grids and thermal model */
% ------------------------------ */

hgrid = (0:10:6000)';   % height grid in feet. MUST be even steps: mixer */
 % in 2007 revision, I'm finding a small grid size is important. How did
 % 500 work before? 
lgrid = 0;                  % thermal strength grid in kts */
lprb  = 1;                  % thermal probability in each mile */
porp  = 0;                  % porpoise factor -- cruse strength/climb */
xmax  = 150;                % course length in n mi */
vwin  = 60;                 % winner's handicapped speed in kts */
vwin = vwin/hcap; 


dx    = 1;                 % step size in n mi */
alpha = 0.65;              % maximum points for landout */
minclimb = 300;            % min altitude for thermaling */
fstp = 1;                  % how many miles back for first step? */
mix = 0;                   % sigma altidude in feet per mile = smoothing  */

                           % Preprogrammed assumptionss */
if assumptions == 1;

    lgrid = [0; (4+smin)];         % thermal strength grid in kts */
    lprb  = [.9; .1];                % thermal probability in each mile */
    % use lgrid = 0; lprb = 1; to check that it reverts to glide to finish */
    porp  = 0;                    % porpoise factor -- cruse strength/climb */
end;

if (assumptions == 2) |(assumptions == 3) |(assumptions == 5)|(assumptions == 7) ;
    lgrid = [0;(1+smin);(2+smin);(4+smin);(6+smin)];
    lprb  = [0;0.2;0.1;0.05;0.02];
    lprb(1) = 1-sum(lprb);
    porp  = 0.5;
    mix = 60; 
end;

if assumptions == 4;
    hgrid = (0:1:9000)';
    lgrid =  [0;   (2+smin);  (4+smin);   (6+smin);(8+smin)];
    lprb  =  [0;0.15;0.1;0.05;0.025];
    lprb(1) = 1-sum(lprb);
    porp  = 0.5;
    vwin = 70;
    xmax = 250;
end;

if (assumptions == 6) ;
    lgrid = [0;(2+smin);(4+smin);(6+smin)];
    lprb  = [0;0.1;0.05;0.02];
    lprb(1) = 1-sum(lprb);
    porp  = 0.5;
    mix = 60; 
end;


% ------------------------------ */
% graphs to assess thermal model */
% ------------------------------ */

if thermgf;
    figure; 
    bar(lgrid,lprb);
    title('Thermal distribution');
    
    figure;
    prbetr = rev(cumsum(rev(lprb)));
    pr1 = 1-prbetr;
    px = (1:1:20)';
    py = 1 - (ones(size(px,1),1)*pr1').^(px*ones(1,size(pr1,1))); 
     % gauss:    py = 1-(pr1')^px;
    plot(px,py);

    title('Chance finding a thermal as good or better than indicated in x miles');
    if (assumptions == 2) ;
        legend('1 kt','2 kts','4 kts','6 kts');
    end;
    if (assumptions == 4) ;
       legend('1 kt','2 kts','4 kts','6 kts','8kts');
    end;
    xlabel('Nautical miles');
    ylabel('Probability');
end;

% ----------------------------- */
% Glider performanc procedures  */
% All procedures use a quadratic*/
% polar approximation           */
% ----------------------------- */

% find coefficient a in s(v) = smin + a/2 (v-vms)^2 to fit s80 */

a = 2*(shi-smin)/(vhi-vms)^2;
ldmax = G(0,0);

% --------------------------------- */
% plot polar and related quantities */
% --------------------------------- */

if polarplt;

    figure;
    subplot(2,2,1); 
    v = (30:1:130)'; 
    plot(v,s(v));
    title('sink rate');
    xlabel('speed in kts');
    axis([35 130 -inf inf]); 
    
    subplot(2,2,2); 
    plot(v,mc(v));
    title('Macready to fly each speed in still air');
    axis([35 130 -inf inf]); 
    
    subplot(2,2,3); 
    plot(v,v./s(v));
    title('glide angle');
    axis([35 130 -inf inf]); 
    
    m = [2 0 -2];
    lam = (0.1:0.1:10)';
    m = ones(size(lam,1),1)*m;
    lam = lam*ones(1,size(m,2)); 
    Ga = G(m,lam);
    
    figure;
    plot(lam,Ga);
    title('Glide as function of MacCready and lift/sink');
    legend('m=+2','m=0','m=-2');
    xlabel('McReady');
    axis([0 10 -100 100]); 
    
    figure; 
    plot(lam,mcinv(lam-m));
    title('Airspeed as function of MacCready and lift/sink');
    legend('m=+2','m=0','m=-2');
    xlabel('MacCready');
end;     

if ~justgrf; % allows you to skip calculations and play with graphs */

% ==================================================================== */
%                      START MACCREADY CALCULATIONS                    */
% ==================================================================== */

% -------------------------------- */
% Find values one mile out         */
% simply fly home, no lift or sink */
% -------------------------------- */

Twin = xmax/vwin; % winner's time -- used to evaluate landouts
                                   % initialize grids */
w = zeros(rows(hgrid),1);          % speed at finish */
wh = zeros(rows(hgrid),1);         % wh (h,x) */
wt = -ones(rows(hgrid),1)/Twin;    % wt (h,x) if no landout */
wt(1) = 0;                         % at zero, you've landed out */
w(1) = alpha*(xmax-fstp*dx)/xmax;
                               % landout rules, assuming winning speedtonow*/
wh(1) = NaN;                  % code for not known yet */
hi = 2;
while hi <= rows(hgrid);
    h = hgrid(hi);
    v = vofg(0,(6000*fstp*dx)/h);     % speed to just exhaust altitude in this lift*/
    if verbose; 
        fprintf('height= %8.3f speed= %8.3f ', h,v) ;
    end;
    if v >= vms;           % can make it home */
        if verbose; 
            fprintf(' can make it home, mc= %8.3f\n', mc(v));
        end;
        w(hi) = xmax/(Twin*(xmax-fstp*dx)/xmax + fstp*dx/v) /vwin;
                           % speed at finish, assuming you've flown at winning speed so far */
        wh(hi) = 1/Twin * 1/ ( mc(v) );
    elseif v == -1;        % can't make it home */
        wt(hi) = 0;
        w(hi) = NaN;
        wh(hi) = NaN;
        if verbose; 
            fprintf('  cant make it home \n');
        end;
        
    else;
        fprintf('error: speed %8.f reported to get home\n',v);
    end;
    hi = hi+1;
end;

if verbose;
    fprintf(' first 10 h  w wh wt before interpolation \n'); 
    fprintf(' %8.8s %8.8s %8.8s %8.8s\n','h','w','wh','wt');
    fprintf(' %8.3f %8.3f %8.3f %8.3f\n', [hgrid(1:10) w(1:10) wh(1:10) wt(1:10)]'); 
end;

% now fill in wh for landout h values by linearly interpolating w */

i = 2;
while isnan(w(i));
    i = i+1;
end;
   % now i is the smallest index of a valid height */
if verbose; 
    fprintf('i index of smallest valid height %8.3f \n', i);
end;

while i <= rows(hgrid);
    testwh = 6000*(w(i)-w(1))/(hgrid(i)-hgrid(1));
                             % Guess at wh in landout region */
    if testwh > wh(i);       % n't want a huge peak in wh */
        break;               % (infinite value of an inch) */
    else;                    % hence use linear w until we reach a */
        i = i+1;             % reasonable wh from the fly home rule */
    end;
end;

wh(1:i-1) = testwh*ones(i-1,1);
w(1:i-1) = w(1)+(hgrid(1:i-1)-hgrid(1))./(hgrid(i)-hgrid(1))*(w(i)-w(1));

if verbose;
    fprintf(' first 10 h  w wh wt after landout,before interpolation \n'); 
    fprintf(' %8.8s %8.8s %8.8s %8.8s\n','h','w','wh','wt');
    fprintf(' %8.3f %8.3f %8.3f %8.3f\n', [hgrid(1:10) w(1:10) wh(1:10) wt(1:10)]'); 
end;

if verbose;
    figure
    subplot(2,2,1)
    plot(hgrid,w);
    title('w one mile out');
    
    subplot(2,2,2)
    plot(hgrid,wh);
    title('wh one mile out');

    subplot(2,2,3); 
    plot(hgrid,wt);
    title('wt one mile out');
    
    subplot(2,2,4); 
    py = -wt./wh;
    px = hgrid*ones(1,cols(py));
    plot(px,py); 
    hold on;
    px = g2(0,py)*fstp*6000;
    plot(px,py)
    legend('calculation','Macready to finish');
    title('MacCready one mile out');
   
end;

% take average over lift values */

whv = wh;
wtv = wt;
lamv = -wtv./whv;
lamfin = lamv*ones(1,fstp); % stores final answer */

% we're done! we have wh and wt one mile out. */

% ----------------  */
% now iterate back  */
% ----------------- */

whl = zeros(rows(hgrid),rows(lgrid));         % wh(h,x,l) */
wtl = zeros(rows(hgrid),rows(lgrid));
whl(1,:) = alpha*ldmax/xmax*ones(1,rows(lgrid)); % landout value */

holdg = hgrid;                    % place to store height at x -dx */
xi = fstp+1;                      % index of distance */
while xi <= xmax+1;

    % mix up h */                % adds subtracts random amount to h -- smooths wh, wt over h */
    whvnew = mixer(whv,hgrid,mix) ;
    wtvnew =  mixer(wtv,hgrid,mix);
    if verbose;
        fprintf('Starting mile %8.3f hit any key\n', xmax-xi) ;% gauss program had  hit any key to continuewait;;
        figure;
        plot(hgrid,[whv whvnew wtv wtvnew]);
        title('old and smoothed wh, wt');
    end;
    whv = whvnew;
    wtv = wtvnew;

    li = 1;                         % look over all lift values */
     while li <= rows(lgrid);
        l = lgrid(li);
        hi = 1;                     % look over all altitudes; hi = h index */
         while hi <= rows(hgrid);
            
            % for each altitude and lift, we have the mcready, so we can work
            %   back and find the altitude at which you leave to get here */            
            holdg(hi) = hgrid(hi) + g2(l*porp,lamv(hi))*dx*6000;
            hi = hi+1;
        end;

        % at this point, whl and wtl are wh and wt, but at the new h.
        % next task is to interpolate this to give wh and wt on the old h grid*/

        hi = 2;   % one always means land out */
         while hi <= rows(hgrid);
            h = hgrid(hi);
            if h <= min(holdg); % you're out */
                whl(hi,li) = alpha*ldmax/xmax;
                wtl(hi,li) = 0;
            else;
                whl(hi,li) = interp(h,holdg,whv);   % interpolates */
                wtl(hi,li) = interp(h,holdg,wtv);
            end;
            hi = hi+1;
        end;

        if verbose;
            figure;
            subplot(2,1,1);
            plot(holdg*[1 1],[whv wtv]);
            hold on
            plot(hgrid*ones(1,2),[whl(:,li) wtl(:,li)]);
            title('1 working back; cruise whl and wtl and interpolation');
            legend('wh before','wt before','wh after','wtafter');

            subplot(2,1,2);
            plot(holdg*(1),(-wtv./whv));
            hold on;
            plot(hgrid*ones(1,1),(-wtl(:,li)./whl(:,li)));
            legend('lam before','lam after');
            title('2 working back; cruise lam and interpolation');            
        end;

        % at this point, we have wht and wlt assuming that you will
        %   glide. However, you have the option to climb. If l-smin > mc
        %   then you climb up to the point that l-smin = mc (mc rises with
        %   altitude) or the top of the lift, whichever comes first. */

         if max(whl(:,li) == 0) ; 
             fprintf('error: whl = 0 whl = \n')
             disp(whl(:,li)');
         end;
         lam = -wtl(:,li)./whl(:,li); % MacCready if you criuse */
         flg = 0;                     % flag for, have you foud point to leave yet */
         if min(lam) < l - smin;  %if there are any altitudes where you climb*/

                                  % find altitude where you stop thermaling */
            indx = rows(hgrid);   % start at top */
             while hgrid(indx) > minclimb; % can't thermal below minclimb */
                if lam(indx) < l - smin;     % if lift too weak, leave here */

                    if flg == 0;             % if this is where you leave */
                        wtall = wtl(indx,li);
                        if indx == rows(hgrid);   % if leave at top, wh jumps*/
                            whall = -wtall/(l - smin);
                        else;
                            whall = -wtall/(l - smin); %let wh jump to soak approx error */
                        end;
                        flg = 1;
                    end;
                   if flg == 1;             % if you will leave higher, Mc work  wn */
                        wtl(indx,li) = wtall;
                        whl(indx,li) = whall;
                    end;
                end;
                indx = indx-1;
            end;
         end;

        li = li+1;                          % go on to next lift */
    end;
    whv = sum( (lprb*ones(1,rows(whl))).*whl',1)';                 % expected values over lift */ % note gauss was (lprb.*whl')
    wtv = sum( (lprb*ones(1,rows(whl))).*wtl',1)';
    lamv = -wtv./whv;
    xi = xi+1;                              % move back one mile */
    lamfin = [lamfin lamv];                   % store final results */

    
end;

if assumptions == 1; 
    save('A1','lamfin'); 
elseif assumptions == 2; 
    save('A2','lamfin'); 
elseif assumptions == 3; 
    save('A3','lamfin'); 
elseif assumptions == 4; 
    save('A4','lamfin'); 
elseif assumptions == 5; 
    save('A5','lamfin'); 
elseif assumptions == 6; 
    save('A6','lamfin'); 
elseif assumptions == 7; 
    save('A7','lamfin'); 
end;


end; % ends justgrf */


% ----------------------------------------------- */
% Graph MacCready rules vs. height, distance to go  */
% ----------------------------------------------- */

% Figures use 2 kts/ms  and 300m/1000ft to keep it simple
if 1;
    load('A2');
    lamfin_discus = lamfin;
    
    figure;
    plot(lamfin_discus((hgrid>=minclimb),[8;16;54;150])/2,hgrid(hgrid>=minclimb)/3,'LineWidth',2);
    
    hold on;
    mcvals = (0:0.1:5)';
    plot(mcvals/2,g2(0,mcvals)*[10 20]*6000/3,'--');
    grid on; 
    set(gca,'Ytick',[0 500 1000 1500 2000])
        
    xlabel('MacCready Wert (m/s)');
    ylabel(['H' char(246) 'he (m)'],'Interpreter','none');
    %legend('15 km','30','100','275');
     
    axis([0 2.5 0 2000]);
    text(2,500,'15 km','Color','b'); 
    text(2,1200,'30 km','Color','g'); 
    text(1.5,1700, '100,','Color','c'); 
    text(1.67,1700,'275 km','Color','r'); 
   % implement as text 
   % _Pmsgstr = '10 nm out\00020 nm out\000100\000150';
   % _pmsgctl =  (4  ~ 1700 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 );
   %             (4  ~ 3800 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 );
   %             (3.2  ~ 5000 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 );
   %             (2.7  ~ 5000 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 );
 % _Pmsgctl =  x(ll) ~ y(ll) ~ ht(.15)  ~angle(o) ~  plot(1)  ~ co(10)~ thick*/
    print -depsc2 fig2.eps
end;


if 1; 
    load('A2');
    lamfin_discus = lamfin;
    load('A5'); 
    lamfin_ka6 = lamfin; 
    load('A6'); 
    lamfin_discus_light = lamfin; 
    load('A7'); 
    lamfin_25 = lamfin; 

    
    figure;
        plot(lamfin_discus(hgrid>=minclimb,150)/2,hgrid(hgrid>=minclimb)/3,'LineWidth',2);
        hold on;
        plot(lamfin_discus_light(hgrid>=minclimb,150)/2,hgrid(hgrid>=minclimb)/3,'--b','LineWidth',2);
        plot(lamfin_ka6(hgrid>=minclimb,150)/2,hgrid(hgrid>=minclimb)/3,'-r','LineWidth',2);
    
        grid on; 
        xlabel('MacCready Wert (m/s)');
        ylabel(['H' char(246) 'he (m)'],'Interpreter','none');
        axis([0 2.5 0 2000]);
        set(gca,'Ytick',[0 500 1000 1500 2000])
        text(1.52,750,'Discus','Color','b'); 
        text(1.2,1200,'Ka6','Color','r'); 
        print -depsc2 fig1.eps

    figure;
        plot(lamfin_discus(hgrid>=minclimb,150)/2,hgrid(hgrid>=minclimb)/3,'LineWidth',2);
        hold on;
        plot(lamfin_discus_light(hgrid>=minclimb,150)/2,hgrid(hgrid>=minclimb)/3,'--b','LineWidth',2);
        plot(lamfin_ka6(hgrid>=minclimb,150)/2,hgrid(hgrid>=minclimb)/3,'-r','LineWidth',2);
        plot(lamfin_25(hgrid>=minclimb,150)/2,hgrid(hgrid>=minclimb)/3,'-m','LineWidth',2);
    
        grid on; 
        xlabel('MacCready Wert (m/s)');
        ylabel(['H' char(246) 'he (m)'],'Interpreter','none');
        axis([0 2.5 0 2000]);
        set(gca,'Ytick',[0 500 1000 1500 2000])
        text(1.2,800,'Discus','Color','b','backgroundcolor','w'); 
        text(1.2,1200,'Ka6','Color','r');
        text(1.5,650,'ASH-25','Color','m'); 

end;