function f_PTplot(T, d, FN, ts)

% plot borehole temperature data at a specified site

%* this script is intended for use with data in a format like 
%  that in PTrecords.mat
%
%  T_*  year  
%       day-of-year 
%       temp1   (this is surface air temperature)
%       temp2 
%       ...
%
%  d_*  depth1  (should be +ve, 2 m)
%       depth2
%       ...
%
%* input variables
%  T:   temperature data set
%  d:   corresponding depths
%  FN:  starting figure number
%  ts:  naming string to identify data site


N=size(T)*[1 0]';   % number of data records
M=size(T)*[0 1]'-2; % number of depths

y1=T(1,1);
yN=T(N,1);

lc=['r- '; 'y- '; 'g- '; 'c- '; 'b- '; 'm- ';
    'r-.'; 'y-.'; 'g-.'; 'c-.'; 'b-.'; 'm-.';
    'r: '; 'y: '; 'g: '; 'c :'; 'b: '; 'm: '];

baddies=find(T==999.9);
T(baddies)=NaN;


%* plot all data at a given depth over entire record
figure(FN)
clf
plot([1:N], T(:,3), 'b.')
hold on
plot([1:N], T(:,4), 'r.')
plot([1:N], T(:,5), 'g.')
plot([1:N], T(:,M), 'm.')
axis tight
ylabel('temperature (C)')
xlabel('cummulative days in record')
title([ts int2str(y1) ' to ' int2str(yN)])
legend('air', num2str(d(4)), num2str(d(5)), num2str(d(M))) 



%* plot every day for one year
tp=find(T(:,1)==y1+1);       % pick second year

figure(FN+1)
clf
for n=tp(1):tp(length(tp))
    plot(T(n,4:M+2), d(2:M))
    hold on
end
xlabel('temperature (C)')
ylabel('depth below the surface (m)')
title([ts int2str(y1+1) ' all days'])




%* plot same days of year at all depths in soil
dtp=30;

figure(FN+2)
clf
for n=y1:yN
    ty=find(T(:,1)==n);
    for t=min(ty):max(ty)
        if t/dtp==round(t/dtp)
            plot(T(t,4:M+2), d(2:M), lc(n-y1+1,:))
            hold on
        end
    end
end
xlabel('temperature (C)')
ylabel('depth below the surface (m)')
title([ts int2str(y1) ' to ' int2str(yN) ' every ' int2str(dtp) ' days'])