average.pro 6.89 KB
; $Id: average.pro,v 1.6 2012/07/06 13:24:01 budnik Exp $
;====================================================================
;--------------- Averaging for 1D & 2D arrays (second dim -> Time) 
;====================================================================

pro AVERAGE, AverStructure, T, Val 
common GraphC, Graph, GraphN, Item

 FillValue = Graph[GraphN].FillValue[Item]
 Yinfo = SIZE(Val) 
 if (Yinfo[0] GT 2) then RETURN
 RetSize = Yinfo[0] eq 2 ? Yinfo[2] : Yinfo[1]

 if (Yinfo[0] eq 2) then index =  Yinfo[1] eq 1 ? 0 : Yinfo[1]/2 ;  
 
 specialCase = 0; 
 
 if (Graph[GraphN].Step GT AverStructure.Sampling*2.0) then begin
           if (finite(FillValue)) then begin

              if (Yinfo[0] eq 1) then num31 = where(Val ne FillValue, Index31) $
              else begin
                 num31 = where(Val[index,0:RetSize-1] ne FillValue, Index31);
 ; artificial truc for different spectra definitions....
                if (index gt 1) then begin                 
                  num31_1 = where(Val[1,0:RetSize-1] ne FillValue, Index31_1);
                  if (Index31_1 ne Index31) then begin                       
                        specialCase = 1;                                                                     
                  endif else begin
                        if (Index31 le 0) then begin
                            numTest =  where(Val[*,0:RetSize-1] ne FillValue, IndexTest);
                            if (IndexTest gt Yinfo[1]) then specialCase = 1;  
                        endif
                  endelse
              endif
              endelse 
           endif else begin
             if (Yinfo[0] eq 1) then num31 = where(finite(Val), Index31) $
             else begin
                num31 = where(finite(Val[index,0:RetSize-1]), Index31); 
; artificial truc for different spectra definitions....
                if (index gt 1) then begin                 
                  num31_1 = where(finite(Val[1,0:RetSize-1]), Index31_1);
                  if (Index31_1 ne Index31) then begin                       
                        specialCase = 1;                                                                     
                  endif else begin
                        if (Index31 le 0) then begin
                            numTest =  where(finite(Val[*,0:RetSize-1]), IndexTest);
                            if (IndexTest gt Yinfo[1]) then specialCase = 1;  
                        endif
                  endelse
              endif
            endelse    
        endelse
   
         if (specialCase) then Index31 = n_elements(T);
 
         if (Index31 le 0) then begin
               AverStructure.Nres = -1L
             return
         endif
  
       if (specialCase) then begin          
            Val = REFORM(Val[*,*],Yinfo[1]*Index31)
       endif else begin 
            T = T[num31]
            Val = Yinfo[0] eq 2 ? REFORM(Val[*,num31],Yinfo[1]*Index31) : Val[num31];
       endelse
   
      MF =  AverStructure.Nres eq -1L ? Val : [*(AverStructure.LastVal),VAL]
      Time = AverStructure.Nres eq -1L ? T : [*(AverStructure.LastTime),T]
      N_Time = N_elements(Time)
      NAve = long(((Time[N_Time-1] < Graph[GraphN].TotalTime) - Time[0]) / Graph[GraphN].Step + 1)> 1L
     
      if Yinfo[0] eq 2 then MF = REFORM(MF, Yinfo[1], N_Time) 

       numStart = where(abs((*AverStructure.Time)-Time[0]) LE Graph[GraphN].Step/2.0, IndexStart)
       if (IndexStart le 0) then begin
            AverStructure.Nres = -1L
            return
       endif
       NStart = numStart[0]
             
      for i = NStart, NStart+NAve-1L  do begin   
         num = where(abs((*AverStructure.Time)[i]-Time) LE Graph[GraphN].Step/2.0, Index)
         if (Index LE 0) AND (Graph[GraphN].Sampling[Item] LT Graph[GraphN].MaxSampling[Item]) then $
              num = where(abs((*AverStructure.Time)[i]-Time) LE Graph[GraphN].MaxSampling[Item]/2.0, Index)
         if (Index GT 0) then begin
          if (Yinfo[0] eq 2) then begin
             if (specialCase) then begin
                 for k = 0, Yinfo[1]-1 do begin
                    if (Index GT 1) then begin
                        temp = 0.0;
                        numK = 0;
                        for kk = 0, Index-1 do begin
                           if  (finite(FillValue)) then begin
                                if (finite(MF[k,num[kk]])) then begin
                                    temp += MF[k,num[kk]];
                                    numK++;
                                endif
                          endif else begin
                                if (MF[k,num[kk]] ne FillValue ) then begin
                                    temp += MF[k,num[kk]];
                                    numK++;
                                endif
                          endelse
                         endfor                          
                        (*AverStructure.Val)[i*Yinfo[1]+k]  = temp/ numK; 
                    endif else (*AverStructure.Val)[i*Yinfo[1]+k] = MF[k,num[0]]
                 endfor
             endif else (*AverStructure.Val)[i*Yinfo[1]:i*Yinfo[1]+Yinfo[1]-1] = Index GT 1 ? total(MF[*,num],2)/Index : MF[*,num[0]] 
          endif else (*AverStructure.Val)[i] = Index GT 1 ? total(MF[num])/Index : MF[num[0]]
         endif  
      endfor  
      AverStructure.Nres = Index GT 0 ? (num[Index-1]+1) < (N_Time-1) : -1L  
      
      if (Graph[GraphN].CurrentSec LT (Graph[GraphN].Sections-1)) AND (AverStructure.Nres ne -1) then begin
            *(AverStructure.LastTime) = Time[AverStructure.NRes:*]
            *(AverStructure.LastVal) = Yinfo[0] eq 2 ? REFORM(MF[*,AverStructure.NRes:N_Time-1], Yinfo[1]*(N_Time-AverStructure.Nres)) :   MF[AverStructure.NRes:N_Time-1] 
      endif
      
      if (Graph[GraphN].CurrentSec eq (Graph[GraphN].Sections-1)) then begin
        if ptr_valid(AverStructure.LastTime) then ptr_free, AverStructure.LastTime
        if ptr_valid(AverStructure.LastVal) then ptr_free, AverStructure.LastVal
      endif
      
   endif else begin ;if (Step GT Sampling*2.0)

     records = N_elements(*AverStructure.Time);
 
     if Yinfo[0] eq 2 then Val = reform(Val, Yinfo[1]*RetSize);
     if (records eq 0) then begin
         (*AverStructure.Time) = T; : [*AverStructure.Time,T]
         (*AverStructure.Val) =  Val;
     endif else begin 
         numBefore = where((*AverStructure.Time)[records-1] GT T, IndexBefore);
         if (IndexBefore eq 0) then begin
                (*AverStructure.Time) =  [*AverStructure.Time,T];
                (*AverStructure.Val) =   [*AverStructure.Val,Val];
         endif else begin
                (*AverStructure.Time) =  [*AverStructure.Time,T[IndexBefore-1:RetSize-1]];
                if (Yinfo[0] eq 1) then  (*AverStructure.Val) =  [*AverStructure.Val,Val[IndexBefore-1:RetSize-1]]
                if (Yinfo[0] eq 2) then  (*AverStructure.Val) =  [*AverStructure.Val,Val[Yinfo[1]*(IndexBefore-1):Yinfo[1]*RetSize-1]]  
         endelse
     endelse
    endelse  

return
end