Home > atmlab > retrieval > qpack2 > qp2_l2.m

qp2_l2

PURPOSE ^

QP2_L2 Standardised L2 output for Qpack2

SYNOPSIS ^

function L2 = qp2_l2( Q, R, xa, O, X, Y, m )

DESCRIPTION ^

 QP2_L2    Standardised L2 output for Qpack2

    The function allows that a L2 structure is generated where the set of
    fields can be controlled through Q. The fields of L2 are determined by the
    settings of the L2 field of the retrieval quantities (such as
    Q.POLYFIT.L2) and Q.L2_EXTRA. The naming and content of the L2 fields
    are described further in the Qpack2 manual.

    Possible choices for L2_EXTRA are:
    ----------: These data are transferred from X, the output of *oem*
         'dx' : Last value of X.dx. See qinfo(@oem,'dx')
       'cost' : Last value of X.cost. See qinfo(@oem,'dx')
          'e' : See qinfo(@oem,'e')
         'eo' : See qinfo(@oem,'eo')
         'es' : See qinfo(@oem,'es')
         'ex' : See qinfo(@oem,'ex')
         'yf' : See qinfo(@oem,'yf')
          'A' : See qinfo(@oem,'A')
          'G' : See qinfo(@oem,'G')
          'J' : See qinfo(@oem,'J')
          'S' : See qinfo(@oem,'S')
         'So' : See qinfo(@oem,'So')
         'Ss' : See qinfo(@oem,'Ss')
    ----------: Other data
       'date' : Date of the measurement. Taken from Y. Gives fields 'year', 
                'hour', 'day', 'hour', 'min' and 'sec'
         'xa' : A priori state for each retrieval quantity.
          'y' : Measured spectrum.
         'bl' : Baseline fit.
     'tnoise' : Thermal noise magnitude. A scalar. Overall mean of Y.TNOISE.
        'ptz' : Includes 'p_grid', 'z_field' and 't_field' from corresponding
                fields of Q. For higher atmospheric dimensionalities 
                'lat_grid' and 'lon_grid* are also included.

    The variables G, J S, So and Ss are stored as complete matrices. The
    matrix A, and the vectors e, eo, es, ex and xa (as well as x and mr)
    are splitted, and are stored separately for each retrieval quantity.

    A field speciesX_vmr0 is added automatically when the retrieval unit is
    rel or logrel.

    It must be ensured that the output of *oem* includes the data to be
    extracted here. This is most easily done by this function, following the
    first format below. Instead of initiating as O = oem;, instead do:
       O = qp2_l2(Q); 
    And O will be initiated in a way that this function can be used for
    producing the L2 data.

 FORMAT   O = qp2_l2( Q )

 OUT   O       OEM setting structure. 
 IN    Q       Qarts setting structure. 

      or

 FORMAT   L2 = qp2_l2( Q, R, xa, O, X, Y, m )

 OUT   L2      Creaed L2 structure.
 IN    Q       Qarts setting structure. 
       R       See *arts_oem_init*.
       O       OEM setting structure. See *oem*.
       xa      A priori state vector.
       X       Output structure from *oem*
       Y       Measurement structure. See qinfo(@qp2_y)
       m       Index of retrieved measurement.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

qp2_l2.m

SOURCE CODE ^

0001 % QP2_L2    Standardised L2 output for Qpack2
0002 %
0003 %    The function allows that a L2 structure is generated where the set of
0004 %    fields can be controlled through Q. The fields of L2 are determined by the
0005 %    settings of the L2 field of the retrieval quantities (such as
0006 %    Q.POLYFIT.L2) and Q.L2_EXTRA. The naming and content of the L2 fields
0007 %    are described further in the Qpack2 manual.
0008 %
0009 %    Possible choices for L2_EXTRA are:
0010 %    ----------: These data are transferred from X, the output of *oem*
0011 %         'dx' : Last value of X.dx. See qinfo(@oem,'dx')
0012 %       'cost' : Last value of X.cost. See qinfo(@oem,'dx')
0013 %          'e' : See qinfo(@oem,'e')
0014 %         'eo' : See qinfo(@oem,'eo')
0015 %         'es' : See qinfo(@oem,'es')
0016 %         'ex' : See qinfo(@oem,'ex')
0017 %         'yf' : See qinfo(@oem,'yf')
0018 %          'A' : See qinfo(@oem,'A')
0019 %          'G' : See qinfo(@oem,'G')
0020 %          'J' : See qinfo(@oem,'J')
0021 %          'S' : See qinfo(@oem,'S')
0022 %         'So' : See qinfo(@oem,'So')
0023 %         'Ss' : See qinfo(@oem,'Ss')
0024 %    ----------: Other data
0025 %       'date' : Date of the measurement. Taken from Y. Gives fields 'year',
0026 %                'hour', 'day', 'hour', 'min' and 'sec'
0027 %         'xa' : A priori state for each retrieval quantity.
0028 %          'y' : Measured spectrum.
0029 %         'bl' : Baseline fit.
0030 %     'tnoise' : Thermal noise magnitude. A scalar. Overall mean of Y.TNOISE.
0031 %        'ptz' : Includes 'p_grid', 'z_field' and 't_field' from corresponding
0032 %                fields of Q. For higher atmospheric dimensionalities
0033 %                'lat_grid' and 'lon_grid* are also included.
0034 %
0035 %    The variables G, J S, So and Ss are stored as complete matrices. The
0036 %    matrix A, and the vectors e, eo, es, ex and xa (as well as x and mr)
0037 %    are splitted, and are stored separately for each retrieval quantity.
0038 %
0039 %    A field speciesX_vmr0 is added automatically when the retrieval unit is
0040 %    rel or logrel.
0041 %
0042 %    It must be ensured that the output of *oem* includes the data to be
0043 %    extracted here. This is most easily done by this function, following the
0044 %    first format below. Instead of initiating as O = oem;, instead do:
0045 %       O = qp2_l2(Q);
0046 %    And O will be initiated in a way that this function can be used for
0047 %    producing the L2 data.
0048 %
0049 % FORMAT   O = qp2_l2( Q )
0050 %
0051 % OUT   O       OEM setting structure.
0052 % IN    Q       Qarts setting structure.
0053 %
0054 %      or
0055 %
0056 % FORMAT   L2 = qp2_l2( Q, R, xa, O, X, Y, m )
0057 %
0058 % OUT   L2      Creaed L2 structure.
0059 % IN    Q       Qarts setting structure.
0060 %       R       See *arts_oem_init*.
0061 %       O       OEM setting structure. See *oem*.
0062 %       xa      A priori state vector.
0063 %       X       Output structure from *oem*
0064 %       Y       Measurement structure. See qinfo(@qp2_y)
0065 %       m       Index of retrieved measurement.
0066 
0067 % 2009-07-01   Created by Patrick Eriksson.
0068 
0069 function L2 = qp2_l2( Q, R, xa, O, X, Y, m )
0070 
0071 if nargin == 1
0072   %
0073   L2 = oem;
0074   %
0075   if ~qarts_isset( Q.L2_EXTRA ) 
0076     return;           % --->
0077   end 
0078   Q.L2_EXTRA = qarts_get( Q.L2_EXTRA );
0079   rqre_datatype( Q.L2_EXTRA, @iscellstr, 'Q.L2_EXTRA' );                    %&%
0080   %
0081   for i = 1 : length( Q.L2_EXTRA )
0082     switch Q.L2_EXTRA{i} 
0083       case {'dx','cost','e','eo','es','ex','yf','A','G','J','S','So','Ss'}
0084         L2.(Q.L2_EXTRA{i})= true;
0085       case 'mresp'
0086         L2.A = true;
0087       case {'date','xa','y','bl','tnoise','ptz'}
0088         %
0089       otherwise
0090         error( sprintf('Unknow out variable: %s', Q.L2_EXTRA{i} ) );
0091     end        
0092   end
0093   %
0094   return;           % --->
0095 end
0096 
0097 %- To handle ase when no putput is selected
0098 %
0099 L2 = [];
0100 
0101 
0102 %- Set up some do-variables to avoid repetition of any(strcmp))
0103 %
0104 do_mresp = any( strcmp( Q.L2_EXTRA, 'mresp' ) );
0105 do_A     = any( strcmp( Q.L2_EXTRA, 'A' ) ) ;
0106 do_xa    = any( strcmp( Q.L2_EXTRA, 'xa' ) );
0107 do_e     = any( strcmp( Q.L2_EXTRA, 'e' ) );
0108 do_eo    = any( strcmp( Q.L2_EXTRA, 'eo' ) );
0109 do_es    = any( strcmp( Q.L2_EXTRA, 'es' ) ) ;   
0110 do_ex    = any( strcmp( Q.L2_EXTRA, 'ex' ) ) ;   
0111 
0112 
0113 % Date and position
0114 %
0115 if any( strcmp( Q.L2_EXTRA, 'date' ) )
0116   L2.year   = Y(m).YEAR;
0117   L2.month  = Y(m).MONTH;
0118   L2.day    = Y(m).DAY;
0119   L2.hour   = Y(m).HOUR;
0120   L2.minute = Y(m).MINUTE;
0121   L2.second = Y(m).SECOND;
0122 end
0123 
0124 
0125 % Retrieval diagnostics
0126 %
0127 L2.converged = X.converged;
0128 %
0129 if ~O.linear
0130   if any( strcmp( Q.L2_EXTRA, 'dx' ) )
0131     L2.dx = X.dx(end);
0132   end
0133 end 
0134 %
0135 if any( strcmp( Q.L2_EXTRA, 'cost' ) )
0136   L2.cost   = X.cost(end);
0137   L2.cost_x = X.cost_x(end);
0138   L2.cost_y = X.cost_y(end);
0139 end
0140 
0141 
0142 % Measurement and fits
0143 %
0144 if any( strcmp( Q.L2_EXTRA, 'y' ) )
0145   if qarts_isset( Y(m).F )
0146     L2.f  = Y(m).F;
0147   elseif Q.SENSOR_DO
0148     L2.f  = qarts_get( Q.SENSOR_RESPONSE_F );
0149   else
0150     L2.f  = qarts_get( Q.F_GRID );    
0151   end
0152   L2.y  = Y(m).Y;
0153 end
0154 if any( strcmp( Q.L2_EXTRA, 'yf' ) )
0155   L2.yf = X.yf;
0156 end
0157 if any( strcmp( Q.L2_EXTRA, 'bl' ) )
0158   L2.bl = R.bl;
0159 end
0160 if any( strcmp( Q.L2_EXTRA, 'tnoise' ) )
0161   L2.tnoise = mean( mean( Y(m).TNOISE ) );
0162 end
0163 
0164 
0165 % PTZ (last forward model fields)
0166 %
0167 if any( strcmp( Q.L2_EXTRA, 'ptz' ) )
0168   L2.p_grid = R.p_grid;
0169   if Q.ATMOSPHERE_DIM >= 2
0170     L2.lat_grid = R.lat_grid;
0171     if Q.ATMOSPHERE_DIM == 3
0172       L2.lon_grid = R.lon_grid;
0173     end
0174   end
0175   L2.t_field = qarts_get( Q.T_FIELD );
0176   L2.z_field = qarts_get( Q.Z_FIELD );
0177 end
0178 
0179 
0180 % Optional matrices
0181 %
0182 if any( strcmp( Q.L2_EXTRA, 'G' ) )
0183   L2.G = X.G;
0184 end
0185 if any( strcmp( Q.L2_EXTRA, 'J' ) )
0186   L2.J = X.J;
0187 end
0188 if any( strcmp( Q.L2_EXTRA, 'S' ) )
0189   L2.S = X.S;
0190 end
0191 if any( strcmp( Q.L2_EXTRA, 'So' ) )
0192   L2.So = X.So;
0193 end
0194 if any( strcmp( Q.L2_EXTRA, 'Se' ) )
0195   L2.Se = X.Se;
0196 end
0197 
0198 
0199 % Start and stop index for each retrieval quantity
0200 I = zeros( length(R.jq), 2 );
0201 for q = 1 : length(R.jq)
0202   I(q,:) = [ R.ji{q}{1}  R.ji{q}{2} ];
0203 end
0204 
0205 
0206 % Measurement response and averaging kernel
0207 %
0208 if do_A & do_mresp 
0209   [mresp,As] = mrespA( X.A, I );
0210 elseif do_A
0211   As = splitA( X.A, I );
0212 elseif do_mresp
0213   mresp = mrespA( X.A, I );
0214 end
0215 
0216 
0217 % Retrieval quantities
0218 %
0219 i_asj    = find([Q.ABS_SPECIES.RETRIEVE]); 
0220 %
0221 for q = 1 : length(R.jq)
0222 
0223   atmfield = false;
0224   ind      = I(q,1) : I(q,2);
0225   siz      = [ length(ind), 1 ];   % Modified below for 2D and 3D variables
0226   
0227   if strcmp( R.jq{q}.maintag, 'Absorption species' )
0228     %--------------------------------------------------------------------------
0229     % It is assumed that abs. species are included first among the ret. quant.
0230     %--------------------------------------------------------------------------
0231     ispecies = i_asj(q);
0232     vname = sprintf( 'Q.ABS_SPECIES(%d)', ispecies );                       %&%
0233     rqre_field( Q.ABS_SPECIES(ispecies), 'L2', vname );                     %&%
0234     out = qarts_get( Q.ABS_SPECIES(ispecies).L2 );
0235     rqre_datatype( out, @isboolean, sprintf('%s.L2',vname) );               %&%
0236     %
0237     if out 
0238       s        = sprintf( 'species%d', q );
0239       atmfield = true;
0240       Astruct  = Q.ABS_SPECIES(ispecies);
0241       %
0242       L2.([s,'_name'])  = R.jq{q}.subtag;
0243     end
0244     
0245   elseif strcmp( R.jq{q}.maintag, 'Atmospheric temperatures' )
0246     %--------------------------------------------------------------------------
0247     rqre_field( Q.T, 'L2', 'Q.T' );                                         %&%
0248     out = qarts_get( Q.T.L2 );
0249     rqre_datatype( out, @isboolean, 'Q.T.L2' );                             %&%
0250     %
0251     if out
0252       s        = 'temperature';
0253       atmfield = true;
0254       Astruct  = Q.T;
0255     end
0256 
0257   elseif strcmp( R.jq{q}.maintag, 'Wind' )
0258     %--------------------------------------------------------------------------
0259     if strcmp( R.jq{q}.subtag, 'u' )
0260       rqre_field( Q.WIND_U, 'L2', 'Q.WIND_U' );                             %&%
0261       out = qarts_get( Q.WIND_U.L2 );
0262       rqre_datatype( out, @isboolean, 'Q.WIND_U.L2' );                      %&%
0263       %
0264       if out
0265         s        = 'wind_u';
0266         atmfield = true;
0267         Astruct  = Q.WIND_U;
0268       end
0269     elseif strcmp( R.jq{q}.subtag, 'v' )
0270       rqre_field( Q.WIND_V, 'L2', 'Q.WIND_V' );                             %&%
0271       out = qarts_get( Q.WIND_V.L2 );
0272       rqre_datatype( out, @isboolean, 'Q.WIND_V.L2' );                      %&%
0273       %
0274       if out
0275         s        = 'wind_v';
0276         atmfield = true;
0277         Astruct  = Q.WIND_V;
0278       end
0279     elseif strcmp( R.jq{q}.subtag, 'w' )
0280       rqre_field( Q.WIND_W, 'L2', 'Q.WIND_W' );                             %&%
0281       out = qarts_get( Q.WIND_W.L2 );
0282       rqre_datatype( out, @isboolean, 'Q.WIND_W.L2' );                      %&%
0283       %
0284       if out
0285         s        = 'wind_w';
0286         atmfield = true;
0287         Astruct  = Q.WIND_W;
0288       end
0289     else                                                                    %&%
0290       error( 'Unknown wind subtag.' );                                      %&%
0291     end
0292     
0293   elseif strcmp( R.jq{q}.maintag, 'Sensor pointing' )
0294     %--------------------------------------------------------------------------
0295     rqre_field( Q.POINTING, 'L2', 'Q.POINTING' );                           %&%
0296     out = qarts_get( Q.POINTING.L2 );
0297     rqre_datatype( out, @isboolean, 'Q.POINTING.L2' );                      %&%
0298     %
0299     if out
0300       s = 'pointing';
0301     end
0302   
0303   elseif strcmp( R.jq{q}.maintag, 'Frequency' )
0304     %--------------------------------------------------------------------------
0305     if strcmp( R.jq{q}.subtag, 'Shift' )
0306       rqre_field( Q.FSHIFTFIT, 'L2', 'Q.FSHIFTFIT' );                       %&%
0307       out = qarts_get( Q.FSHIFTFIT.L2 );
0308       rqre_datatype( out, @isboolean, 'Q.FSHIFTFIT.L2' );                   %&%
0309       if out
0310         s = 'fshift';
0311       end
0312     elseif strcmp( R.jq{q}.subtag, 'Stretch' )
0313       rqre_field( Q.FSTRETCHFIT, 'L2', 'Q.FSTRETCHFIT' );                   %&%
0314       out = qarts_get( Q.FSTRETCHFIT.L2 );
0315       rqre_datatype( out, @isboolean, 'Q.FSTRETCHFIT.L2' );                 %&%
0316       if out
0317         s = 'fstretch';
0318       end
0319     else                                                                    %&%
0320       error( 'Unknown frequency subtag' );                                  %&%
0321     end
0322     
0323   elseif strcmp( R.jq{q}.maintag, 'Polynomial baseline fit' )
0324     %--------------------------------------------------------------------------
0325     rqre_field( Q.POLYFIT, 'L2', 'Q.POLYFIT' );                             %&%
0326     out = qarts_get( Q.POLYFIT.L2 );
0327     rqre_datatype( out, @isboolean, 'Q.POLYFIT.L2' );                       %&%
0328     %
0329     if out
0330       s = sprintf( 'polyfit%d', sscanf(R.jq{q}.subtag(end+[-1:0]),'%d') );
0331     end
0332 
0333   elseif strcmp( R.jq{q}.maintag, 'Sinusoidal baseline fit' )
0334     %--------------------------------------------------------------------------
0335     rqre_field( Q.SINEFIT, 'L2', 'Q.SINEFIT' );                             %&%
0336     out = qarts_get( Q.SINEFIT.L2 );
0337     rqre_datatype( out, @isboolean, 'Q.SINEFIT.L2' );                       %&%
0338     %
0339     if out
0340       s = sprintf( 'sinefit%d', sscanf(R.jq{q}.subtag(end+[-1:0]),'%d')+1 );
0341     end
0342     
0343   else
0344     error( sprintf( 'Unknown retrieval quantity: %s', R.jq{q}.maintag ) );
0345   end
0346 
0347   if out
0348 
0349     %--- Special stuff for atmospheric fields
0350     if atmfield
0351 
0352       % Output follows strictly retrieval range (empty L2_RANGE)
0353       if ~isfield( Astruct, 'L2_RANGE' )  
0354         %
0355         % Fill special L2 fields
0356         L2.([s,'_p'])     = Astruct.GRIDS{1};
0357         %
0358         if Q.ATMOSPHERE_DIM >= 2
0359           L2.([s,'_lat']) = Astruct.GRIDS{2};
0360           siz(2)          = length( Astruct.GRIDS{2} );
0361           %
0362           if Q.ATMOSPHERE_DIM == 3
0363             L2.([s,'_lon']) = Astruct.GRIDS{3};
0364             siz(3)          = length( Astruct.GRIDS{3} );
0365           end
0366         end      
0367       
0368       % Modify range?
0369       else
0370         %- Check that L2_RANGE is OK
0371         range = qarts_get( Astruct.L2_RANGE );
0372         rqre_datatype( range, @isvector, sprintf('L2 range for %s',s) );    %&%
0373         if length(range) < Q.ATMOSPHERE_DIM*2                               %&%
0374           error( sprintf('L2 range for %s is too short.',s) );              %&%
0375         end                                                                 %&%
0376         if range(1) >= range(2)                                             %&%
0377           error( sprintf( ...                                               %&%
0378              'Pressure limits in L2 range (for %s) in wrong order? .',s) ); %&%
0379         end                                                                 %&%
0380 
0381         % Pressure dimension
0382         [ngrid,ip]    = r2g_local( Astruct.GRIDS{1}, R.p_grid, range(1:2) );
0383         siz(1)        = length(ngrid);    
0384         L2.([s,'_p']) = ngrid;
0385         %
0386         % Latitude dimension
0387         if Q.ATMOSPHERE_DIM == 1
0388           ind              = ind( ip );
0389           Astruct.GRIDS{1} = ngrid;     % Used to get VMR0 below
0390         elseif Q.ATMOSPHERE_DIM >= 2
0391           [ngrid,ilat] = r2g_local( Astruct.GRIDS{2}, R.lat_grid, range(3:4) );
0392           siz(2)       = length(ngrid);    
0393           L2.([s,'_lat']) = ngrid;
0394           %
0395           if Q.ATMOSPHERE_DIM == 2
0396             ind = ind( matvec( repmat( ip, 1, length(ilat) ) + ...
0397                                repmat( (ilat'-1)*length(Astruct.GRIDS{1}), ...
0398                                        length(ip), 1 ) ) );
0399             Astruct.GRIDS{2} = ngrid;
0400           else
0401             [ngrid,ilon] = r2g_local( Astruct.GRIDS{3}, R.lon_grid, ...
0402                                                                   range(5:6) );
0403             siz(3)       = length(ngrid);    
0404             L2.([s,'_lon']) = ngrid;
0405             ind = ind( matvec( repmat( ip, [1 length(ilat) length(ilon)] ) + ...
0406                                repmat( (ilat'-1)*length(Astruct.GRIDS{1}), ...
0407                                        [length(ip) 1 length(ilon)] ) + ...
0408                                repmat( (reshape(ilon,[1 1 length(ilon)])-1)*...
0409                                           length(Astruct.GRIDS{1}) * ...
0410                                           length(Astruct.GRIDS{2}), ...
0411                                        [1 length(ilat) length(ilon)] ) ) );
0412             Astruct.GRIDS{3} = ngrid;
0413           end
0414         end
0415         clear ngrid ip ilat ilon
0416       end
0417     end
0418   
0419     %--- Fill L2
0420     %
0421     L2.([s,'_x'])  = reshape( X.x(ind), siz );
0422     %
0423     if any( strcmp( Q.L2_EXTRA, 'xa' ) )
0424       L2.([s,'_xa']) = reshape( xa(ind), siz );
0425     end
0426     %
0427     if any( strcmp( Q.L2_EXTRA, 'e' ) )
0428       L2.([s,'_e'])  = reshape( X.e(ind), siz );
0429     end
0430     if any( strcmp( Q.L2_EXTRA, 'eo' ) )
0431       L2.([s,'_eo']) = reshape( X.eo(ind), siz );
0432     end
0433     if any( strcmp( Q.L2_EXTRA, 'es' ) )    
0434       L2.([s,'_es']) = reshape( X.es(ind), siz );  
0435     end
0436     if any( strcmp( Q.L2_EXTRA, 'ex' ) )    
0437       L2.([s,'_ex']) = reshape( X.ex(ind), siz );  
0438     end
0439     %
0440     if any( strcmp( Q.L2_EXTRA, 'mresp' ) )
0441       L2.([s,'_mr']) = reshape( mresp(ind), siz );
0442     end
0443     if any( strcmp( Q.L2_EXTRA, 'A' ) )
0444       L2.([s,'_A']) = As{q};
0445     end
0446     
0447     % If ABS_SPECIES + rel or logrel, save also vmr0
0448     if strcmp( R.jq{q}.maintag, 'Absorption species' )
0449       if any( strcmp( R.jq{q}.mode, { 'rel', 'logrel' } ) )
0450         vmr0 = arts_regrid( Q.ATMOSPHERE_DIM, Q, getdims( ...
0451                      R.vmr_field0(ispecies,:,:,:), 2:Q.ATMOSPHERE_DIM+1 ), ...
0452                             Astruct.GRIDS );
0453         L2.([s,'_vmr0']) = reshape( vmr0, siz );  
0454       end    
0455     end
0456   end
0457 end
0458 
0459 
0460 
0461 
0462 
0463 function [newgrid,ii] = r2g_local( retgrid, fmgrid, range )
0464   %
0465   np    = length( retgrid );
0466   %
0467   if range(1) >= retgrid(end)  &  range(2) <= retgrid(1)
0468     ii      = find( retgrid >= range(1)  &  retgrid <= range(2) );
0469     newgrid = retgrid(ii);     
0470   elseif range(1) < retgrid(end)  &  range(2) > retgrid(1)
0471     newgrid = [ min([range(2),fmgrid(1)]); 
0472                 retgrid; 
0473                 max([range(1),fmgrid(end)]) ];     
0474     ii      = [ 1 1:np np ]';
0475   elseif range(1) < retgrid(end) 
0476     ii      = find( retgrid <= range(2) );
0477     newgrid = [ retgrid(ii); max([range(1),fmgrid(end)]) ];     
0478     ii      = [ ii; np ];
0479   else
0480     ii      = find( retgrid >= range(1) );
0481     newgrid = [ min([range(2),fmgrid(1)]); retgrid(ii) ];     
0482     ii      = [ 1; ii ];
0483   end
0484 return

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005