function [varargout] = water_temp(varargin)
%This script will measure the temperature of water molecules
%Inputs are dump file, initial timestep and lammps log file
%Output is average temperature of water

%Read temp data from lammps
done = 0;
lammpsTemp = zeros(1,99);
k = 1;
dump = fopen(varargin{3},'r');
while done == 0
    id = fgetl(dump);
    if (strcmpi(id,'START'))
        id = fgetl(dump);
        id = fgetl(dump);
        while k < 100
            data = sscanf(id, '%f', [1 7]);
            tempL = data(2);
            lammpsTemp(k) = tempL;
            id = fgetl(dump);
            k = k+1;
            disp(k);
        end
        done = 1;
    end
end
fclose(dump);

try
    dump = fopen(varargin{1},'r');
catch
    error('Dumpfile not found!');
end

%Read timestep from user input, initialise variables
userstep = varargin{2};
oxygen = 1;
hydrogen = 2;
B = 1.38064852e-23;         %boltzmann's constant
A = 6.0221413e+23;          %avogadro's constant

steps = 99;
temp = zeros(1,steps);
done = 0;
n = 0; 
m = 1;
z = 0;
q = 0;
bound = 100000;

%compute conversions
mO = 16.0/A;
mH = 1.0/A;
kO = (mO*(10^7))/(3*B);
kH = (mH*(10^7))/(3*B);

%Evaluate Temperature from atom velocities 
while done == 0
    %Take a new line in dump file
    id = fgetl(dump);
    
    %Look for a timestep
    %Find timestep
    if (strcmpi(id,'ITEM: TIMESTEP'))
        id = fgetl(dump);
        dumpstep = str2double(id);
        disp(dumpstep);
    end
    
    %check if timestep matches user selection
    if dumpstep == userstep
        %Skip 7 lines to get to atom data
        id = fgetl(dump);
        id = fgetl(dump);
        id = fgetl(dump);
        id = fgetl(dump);
        id = fgetl(dump);
        id = fgetl(dump);
        id = fgetl(dump); %sitting on 'ITEM:ATOMS...'
        
        %Start loop until end of timestep of interest
            while 1
                %Step line in dataset
                id = fgetl(dump);
                
                %At end of timestep of interest break loop
                if strcmpi(id,'ITEM: TIMESTEP')
                    %compute temp for previous timestep
                    
                    % work out average KE = 0.5mv^2/n
                    temp(m) = (0.5*sum(countO)*mO + 0.5*sum(countH)*mH)/(n+q);
                    
                    %find temp from (2/3 KE)/B *units conversion
                    temp(m) = (temp(m)*2*10^7)/(3*B);
                    
                    id = fgetl(dump);
                    dumpstep = str2double(id);
                    disp(dumpstep);
                    if dumpstep < bound
                        m = m+1;
                        n = 0;
                        q = 0;
                        countO = zeros(1,2516);
                        c = 1;
                        while c < 9
                            id = fgetl(dump);
                            c = c+1;
                        end
                        %disp(id);
                    else
                        break
                    end
                end
                
                %convert data into vector
                data = sscanf(id, '%f', [1 9]);
                
                %read atom type
                type = data(3);
                
                %check atom type
                if type == oxygen
                    %read atom velocity
                    z = data(6);
                    vx = data(7);
                    vy = data(8);
                    vz = data(9);
                    if z < 200
                        %evaluate v^2 for atom
                        v2O = (vx*vx)+(vy*vy)+(vz*vz);
                        
                        %write data to vector
                        n = n+1;
                        countO(n) = v2O;
                    end
                end
                
                if type == hydrogen
                    z = data(6);
                    vx = data(7);
                    vy = data(8);
                    vz = data(9);
                    if z < 200
                        %evaluate v^2 for atom
                        v2H = (vx*vx)+(vy*vy)+(vz*vz);
                        %write data to vector
                        q = q+1;
                        countH(q) = v2H;
                    end
                end
                  
            end
            %finish programme
            done = 1;
    end
end

t = linspace(2,200,steps);
plot(t,lammpsTemp,t,temp)
xlabel('Time (ps)')
ylabel('Temperature (K)')
title('Lammps log Temperature compared to Matlab Temperature')
legend('Lammps','Matlab')

fclose(dump);