reverse-shooting

Matlab scripts for reverse shooting
Log | Files | Refs | README

interplin5.m (1673B)


      1 function newy=interplin5(y,x,newx,useclosest)
      2 
      3 % interplin5   5/24/04 -- extrapolates for outside the range
      4 %     function newy=interplin(y,x,newx)
      5 %
      6 %   Simple linear interpolation, assuming that x is monotone
      7 %
      8 %     y = Variable to be interpolated  -- Nx1
      9 %     x = Values we know, corresponding to existing y -- Nx1
     10 %   newx = x values for which we need y values (can include x)
     11 %                    1960 10
     12 %                    1965 15
     13 %                    1967 14
     14 %
     15 %    newy=interplin(y,x,[1962; 1966]) will return [12; 14.5]
     16 %
     17 %   5/24/04:  If you ask for the 1959 value, it will now extrapolate from
     18 %   linearly from the two closest observations.
     19 %
     20 % 7/7/04: If useclosest==1, then will use closest value for extrapolating
     21 %         rather than extrapolating linearly. (Still interpolates in interior).
     22 
     23 if exist('useclosest')==0; useclosest=0; end;  % default is linear interp.
     24 
     25 newy=zeros(length(newx),1);
     26 for i=1:length(newx);
     27    lessthan=find(x<=newx(i));
     28    greatthan=find(x>=newx(i));
     29    jlow=max(lessthan);
     30    jhigh=min(greatthan);
     31    
     32    if useclosest==0;
     33      if isempty(lessthan);  % Extrapolate backwards
     34        jlow=1; jhigh=2;
     35      end;
     36      if isempty(greatthan); % Extrapolate forwards
     37        jhigh=length(x); jlow=jhigh-1;
     38      end;
     39    else;   % use the closest observation at ends rather than linear interp.
     40      if isempty(lessthan);
     41        jlow=1; jhigh=1;
     42      end;
     43      if isempty(greatthan); % Extrapolate forwards
     44        jhigh=length(x); jlow=jhigh;
     45      end;
     46    end;
     47    
     48    if jlow==jhigh; 
     49       newy(i)=y(jlow); 
     50    else;
     51       dx=x(jhigh)-x(jlow);
     52       dy=y(jhigh)-y(jlow);
     53       newy(i)=y(jlow)+dy/dx*(newx(i)-x(jlow));
     54    end;
     55 end;
     56