function hm_temp_movies
%
%   hm_temp_movies
%
%   Function to create a movie of temperature passing by the horizontal
%   array at site B (MBIWE) for each internal wave event.
%   INCLUDES FSIS

%Soupy Alexander, 3/31/04

%This Matlab m-file was used to create portions of U.S. Geological Survey 
%Data Series 85. Although this program has been used by the USGS, no 
%warranty, expressed or implied, is made by the USGS or the United States 
%Government as to the accuracy and functioning of the program and related 
%program material nor shall the fact of distribution constitute any such 
%warranty, and no responsibility is assumed by the USGS in connection 
%therewith.

%Citation:  Butman, Bradford, Alexander, P. Soupy, Anderson, S.P., 
%Lightsom, F.L., Scotti, Alberto, and Beardsley, R.C., 2004, The 
%Massachusetts Bay Internal Wave Experiment, August 1998:  Data Report:  
%U.S. Geological Survey Data Series 85, 1 DVD-ROM.

dataDir = 'C:\SCHTUFF\MASS_BAY\MBAY_IWAVE\REPORT\DATAFILES';

p1Files = {'sbe0035.epic'; 'sbe0052.epic'; 'sbe0045.epic'; ...
    'sbe0048.epic'; 'sbe0044.epic'; 'tpod8093n.epic'};
p1Depths = [20 25 30 35 40 45];
p1Loc = 20;
p2Files = {'sbe0040.epic'; 'sbe0007.epic'; 'sbe0054.epic'; ...
    'sbe0046.epic'; 'sbe0050.epic'; 'tpod8094n.epic'};
p2Depths = [20 25 30 35 40 45];
p2Loc = 50;
p3Files = {' '; 'mcat0683.epic'; 'mcat0686.epic'; 'mcat0685.epic'; ' '; ...
    'mcat0684.epic'};
p3Depths = [NaN 25 30 35 NaN 45];
p3Loc = 80;
p4Files = {'sbe0049.epic'; 'sbe0051.epic'; ' '; 'sbe0039.epic'; ...
    'sbe0101.epic'; 'tpod8091n.epic'};
p4Depths = [20 25 NaN 35 40 45];
p4Loc = 110;
p5Files = {'sbe0041.epic'; 'sbe0038.epic'; 'sbe0047.epic'; ...
    'sbe0103.epic'; ' '; 'tpod8090n.epic'};
p5Depths = [20 25 30 35 NaN 45];
p5Loc = 140;
p6Files = {'fsi1467-a.nc'};
p6Depths = [20];
p6Loc = 25;
p7Files = {'fsi1468-a.nc'};
p7Depths = [20];
p7Loc = 80;
p8Files = {'fsi1469-a.nc'};
p8Depths = [20];
p8Loc = 135;

siteDepths = [p1Depths(:); p2Depths(:); p3Depths(:); ...
    p4Depths(:); p5Depths(:)];
siteLocs = [repmat(p1Loc,6,1); repmat(p2Loc,6,1); repmat(p3Loc,6,1); ...
    repmat(p4Loc,6,1); repmat(p5Loc,6,1)];

depthVec = [20:5:45];
distVec = [20:5:140];

depthOrig = [20:5:45];
distOrig = [20 50 80 110 140];

timeVec = julian(1998, 8, 6, 0):1/(24*60):julian(1998, 9, 3, 0);

tempOrig = repmat(NaN, [length(depthOrig) length(distOrig) length(timeVec)]);

for indexLoc = 1:8;
    eval(['theFiles = p' num2str(indexLoc) 'Files;']);
    eval(['theDepths = p' num2str(indexLoc) 'Depths;']);

    for indexDepth = 1:length(theDepths)
        if isnan(theDepths(indexDepth))
            continue
        end

        ncID = netcdf(fullfile(dataDir, theFiles{indexDepth}), 'nowrite');
        theTime = singlejd(ncID{'time'}(:), ncID{'time2'}(:));
        if ismember(indexLoc, [6 7 8])
            theTime = theTime - 45/(24*60*60);
        end

        if ismember(indexLoc, [1 2 3 4 5])

            theTemp = ncID{'temp'}(:);
            if isempty(theTemp)
                theTemp = ncID{'temperature'}(:);
            end
            theTemp(find(theTemp > 1e3 | theTemp < -1e3)) = NaN;
            tempInterp = interp1(theTime, theTemp, timeVec);
            tempOrig(indexDepth, indexLoc, :) = tempInterp;

        else 
            thePost = find(indexLoc == [6 7 8]);
            velW_this = ncID{'vertv'}(:);
            velW_this(find(velW_this > 1e3 | velW_this < -1e3)) = NaN;
            velU_this = ncID{'east'}(:);
            velU_this(find(velU_this > 1e3 | velU_this < -1e3)) = NaN;
            velW_this = interp1(theTime, velW_this, timeVec);
            velU_this = interp1(theTime, velU_this, timeVec);
            velW{thePost} = -1*velW_this;   %Correct for the plot in axis ij
            velU{thePost} = velU_this;
        end
    end
end

velDepths = [20 20 20];
velLocs = [25 80 135];

[distOrigB, depthOrigB] = meshgrid(distOrig, depthOrig);
[distVecB, depthVecB] = meshgrid(distVec, depthVec);

eventFile = 'C:\SCHTUFF\MASS_BAY\MBAY_IWAVE\REPORT\OTHER\arrival_times.xls';
[theNums, theText] = xlsread(eventFile);

theDates = theText(3:end-2,2);
eventTimes = datenum2julian(datenum(theDates) + theNums(1:end-2,3));
eventData = cell(0,0);
eventTimeVec = cell(0,0);
eventAngle = cell(0,0);
eventSpeed = cell(0,0);

for indexEvent = 2:length(eventTimes);
    fprintf(['On event number ' num2str(indexEvent) ...
        ' of ' num2str(length(eventTimes)) '.\n']);
    eventStart = eventTimes(indexEvent);
    eventEnd = eventStart + 2/24;   %2 hour events
    timeVecEvent = eventStart:1/(24*60):eventEnd;    %1 minute data
    
    inEvent = find(timeVec >= eventStart & timeVec <= eventEnd);
    
    thisEventOrig = tempOrig(:,:,inEvent);
    
    thisEventNew = repmat(NaN, [length(depthVec) ...
        length(distVec) length(inEvent)]);
    
    for timeStep = 1:length(inEvent);
        thisTime = squeeze(thisEventOrig(:,:,timeStep));
        thisTimeOld = thisTime;
        [badRows, badCols] = find(isnan(thisTime));
        for indexBad = 1:length(badRows)
            rowMin = max([1 badRows(indexBad)-1]);
            rowMax = min([size(thisTime, 1) badRows(indexBad)+1]);
            colMin = max([1 badCols(indexBad) - 1]);
            colMax = min([size(thisTime, 2) badCols(indexBad)+1]);
            thisTime(badRows(indexBad), badCols(indexBad)) = ...
                nanmean(nanmean(thisTime(rowMin:rowMax,colMin:colMax)));
        end
    warning off MATLAB:griddata:DuplicateDataPoints
        thisTimeNew = griddata(distOrigB, depthOrigB, ...
            thisTime, distVecB, depthVecB);
        thisEventNew(:,:,timeStep) = thisTimeNew;
    end
    
    eventData{indexEvent} = thisEventNew;
    eventTimeVec{indexEvent} = timeVec(inEvent);
    
    for currSpot = 1:3;
        hitW = velW{currSpot};
        hitW = hitW(inEvent);
        hitU = velU{currSpot};
        hitU = hitU(inEvent);
        [hitAngle, hitRad] = cart2pol(hitU, hitW);
        hitAngle = hitAngle.*180/pi;
        currData(currSpot).angles = hitAngle;
        currData(currSpot).radius = hitRad;
    end
    eventCurr{indexEvent} = currData;
end

%Do the actual movie creation
for indexEvent = 2:length(eventData)
    
    movieFile = ...
        avifile(['C:\SCHTUFF\MASS_BAY\MBAY_IWAVE\REPORT\MOVIES\hmtemp_' ...
        num2str(indexEvent) '.avi'], 'compression', 'none', 'fps', 2);
    
    thisEvent = eventData{indexEvent};
    thisTimeVec = eventTimeVec{indexEvent};
    thisCurr = eventCurr{indexEvent};
    
    figure, hold on
    %set(gcf, 'position', [39 558 1028 286])
    %set(gcf, 'position', [39 558 1028/2 286/2])
    set(gcf, 'position', [39 558 1028*2/3 286*2/3])
    axis([-144 -16 16 46])
    axis ij
    box('on')
    xlabel('Distance From East Buoy (m)')
    ylabel('Depth (m)')
    theXs = get(gca, 'xticklabels');
    theXs(:,1) = [];
    set(gca, 'xticklabels', theXs);
    set(gca, 'dataaspectratiomode', 'man')
    set(gca, 'dataaspectratio', [1 1 1])
    plot(-1*siteLocs, siteDepths, 'k*')
    set(gca, 'clipping', 'off')
    
    velScale = 2;
    l = arrowsafe(-130, 40, 10/velScale, 0);
    set(l, 'color', 'k', 'linewidth', 2)
    t = text(-129.5 + 5/velScale, 42, '10 cm/s');
    set(t, 'horiz', 'center')
    
    for indexTime = 1:size(thisEvent,3)
        thisTime = thisEvent(:,:,indexTime);
        p = pcolor(-1.*distVecB, depthVecB, thisTime);
        caxis([5 12])
        shading interp
        colorbar
        
        theTime = julian2datenum(thisTimeVec(indexTime));
        title(['HM Temperature, Event ' num2str(indexEvent) ', ' ...
            datestr(theTime, 0)]);
        
        for indexCurr = 1:3
            a(indexCurr) = arrowsafe(-1*velLocs(indexCurr), ...
                velDepths(indexCurr), ...
                (thisCurr(indexCurr).radius(indexTime))/velScale, ...
                thisCurr(indexCurr).angles(indexTime));
            set(a(indexCurr), 'color', 'k', 'linewidth', 2, 'clipping', 'off')
        end
        
        
        F = getframe(gcf);
        movieFile = addframe(movieFile, F);
        delete(p)
        delete(a(1))
        delete(a(2))
        delete(a(3))
    end
    movieFile = close(movieFile);
    clear movieFile
    close
end