%  col_prof.m
%  
%  plots a color profile of the temperature or a constituent
%  run get_xtsr first
%  Written by Jim Bowen 7/1/97

%  define the constituent names

cname=['Tracer     ';'Tot Sus Sol';'Coliform   ';'Salinity   ';'Labile DOM ';...
       'Refract DOM';'Algae      ';'Labile POM ';'Ortho P    ';'Ammonia    ';...
       'NO3 + NO2  ';'Diss Oxygen';'Sediment   ';'Totl CO2   ';'Alkalinity ';...
       'pH         ';'CO2 + H2CO3';'Bicarbonate';'Carbonate  ';'Iron       ';...
       'Carbona BOD';'Temperature'];
     
% clear movie frames

clear M

%  set the mpgwrite options, see mpgwrite.doc for an explanation

options = [1, 0, 1, 0, 20, 16, 20, 25];


% compute downstream distances from mouth

for i=1:nseg
  dndis(i)=57.3-0.001*xdist(i);
end


go_on='n';
while go_on=='n'
  disp(' Constituent numbers available')
  disp(kprn)
  cp = input(' Choose a constituent:  ');
% find the constituent number


  cn=1;
  while (kprn(cn) ~= cp)
     cn=cn+1;
  end

%  check if correct constituent was chosen

  disp(['Chosen constituent is ',cname(cp,:)])
  go_on=input(' OK to go on? (y or n):  ','s');
  
end

%  choose output option

outop='x';

while outop~='m' & outop~='i' & outop~='b'
  outop=input('Do you want a movie, images, or both? (m, i, or b): ','s');
end
if outop =='m' | outop=='b'
  disp('Creating movies')
end
if outop =='i' | outop=='b'
  disp('Creating images')
end



%  find the x and z data limits

  xmin=min(dndis); xmax= max(dndis);
  
  zmin=min(min(min(c(:,:,:,cn))))
  zmax=max(max(max(c(:,:,:,cn))))
  
  zc=input(' Do you want an alternate value for zmin? (y or n):  ','s');
  if zc == 'y'
  	zmin=input('Enter alternate value for zmin:  ');
  end
  zc=input(' Do you want an alternate value for zmax? (y or n):  ','s');
  if zc == 'y'
  	zmax=input('Enter alternate value for zmax:  ');
  end

%  set the color map for this constituent

map = jet;
if cp == 12
	map=flipud(jet);
end

%   open figure and set size
figure
set(gcf,'Position',[580,520,560,330])

%  define .mpeg filename

mpeg_file=[infile(1:6),'_',int2str(cp),'.mpg'];


% create a matrix for the movie frames

if outop=='m' | outop=='b'
  M=zeros(length(getframe(gcf)),num_times);
end

% initialize data min and max

zmin_data = 999999.9;
zmax_data = -999999.9;

% set the elevation min and max

ymin=-13.0;
ymax=1.0;

%  start the time counter

for jp=1:num_times

% do only the first and last for outop=='i'

if jp==1 | jp==num_times | outop~='i'

  strtim=num2str(jday(jp));

%   start the data counter

  j=0;
  clear x y z tri trig

%  load the data to be plotted into the  x, y, and z arrays
% do each segment

  for i=1:nseg

%  do the surface layers

    j=j+1;
    x(j)=dndis(i);
    y(j)=el(i,jp);
    z(j)=c(kt(jp),i,jp,cn);
    kx(j)=i;

% and the others

    for k=kt(jp)+1:kb(i)
      j=j+1;
      x(j)=dndis(i);
      y(j)=elb(k);
      z(j)=c(k,i,jp,cn);
      kx(j)=i;
    end
  end

%  do dealaunay triangulation

  tri = delaunay(x,y);

%  set the rest of the limits

 
%  ymin= min(y); ymax= max(y);
  xtext=xmin+0.1*(xmax-xmin);
  ytextt=ymin+0.15*(ymax-ymin);
  ytextb=ymin+0.1*(ymax-ymin);
  zmin_data=min(zmin_data,min(z));
  zmax_data=max(zmax_data,max(z));
  
% count bad triangles

  A= size(tri);
  num_bad=0;
  num_tri=A(1);
  for nt=1:num_tri
    if abs(kx(tri(nt,1))-kx(tri(nt,2)))>1 
       num_bad=num_bad+1;
       bad_tri(num_bad)=nt;
    elseif abs(kx(tri(nt,2))-kx(tri(nt,3)))>1 
       num_bad=num_bad+1;
       bad_tri(num_bad)=nt;
    elseif abs(kx(tri(nt,1))-kx(tri(nt,3)))>1 
       num_bad=num_bad+1;
       bad_tri(num_bad)=nt;
    end
  end

% eliminate bad triangles

  kbad=1;
  kg=0;
  for nt=1:num_tri
      while bad_tri(kbad) < nt
	kbad=kbad+1;
      end
      if nt ~= bad_tri(kbad)
	kg=kg+1;
	trig(kg,:)=tri(nt,:);
      end
  end

%  make the contour plot

  colormap(map);
  trisurf(trig,x,y,z);
  view(0,90);
  axis([xmin xmax ymin ymax]);
  caxis([zmin zmax]);
  set(gca,'Position',[.10,.15,.82,.75])
  text(xtext,ytextt,cname(cp,:))
  text(xtext,ytextb,['Neuse River Profile, jday=',strtim,', run= ',infile])

%  add the location legend and label

 
 xl1=36; yl1=-8.; yli=0.65 ;
  text(xl1,yl1,'SF = Streets Ferry')
 yl1=yl1-yli;
  text(xl1,yl1,'NB = New Bern')
 yl1=yl1-yli;
  text(xl1,yl1,'HS = Hampton Shoal')
 yl1=yl1-yli;
  text(xl1,yl1,'CP = Cherry Point')
 yl1=yl1-yli;
  text(xl1,yl1,'OR = Oriental')
 
 yl=1.3;
  text(dndis(1),yl,'SF','FontWeight','Bold')
  text(dndis(23),yl,'NB','FontWeight','Bold')
  text(dndis(29),yl,'HS','FontWeight','Bold')
  text(dndis(32),yl,'CP','FontWeight','Bold')
  text(dndis(35),yl,'OR','FontWeight','Bold')
  xlabel('Distance downstream from Streets Ferry (km)')
  ylabel('Elevation (m)')
  shading interp
  colorbar

%  save the figure to a movie frame

if outop=='m' | outop=='b'
  M(:,jp)=getframe(gcf);
end

%  create a digital image of the figure

if outop~='m' & (jp == 1 | jp==num_times )
  [XX,map1] = capture(gcf);
  imwrite(XX,map1,'image.pcx')
  if jp==1
    !alchemy -g image.pcx f.gif
  else
    !alchemy -g image.pcx l.gif
  end
  !rm image.pcx
    
end
 
%  end the if statement and the time loop

end
end

disp('Minimum of plotted data')
disp(zmin_data)
disp('Maximum of plotted data')
disp(zmax_data)

%  create an .mpeg movie and save to disk

if outop=='m' | outop=='b'
  disp([' Creating .mpeg movie - Saving file ',mpeg_file])
  mpgwrite(M,map,mpeg_file,options);
end
disp('Done')