clear all
close all

%filename=['darwin_transect_surface_may.nc'];
%filename=['darwin_transect_surface_allmonths.nc'];
%filename=['darwin_transect_0to250m_may.nc'];
 filename=['darwin_transect_0to250m_allmonths.nc'];

lat=ncread(filename,'lat');
depth=ncread(filename,'depth');
month=ncread(filename,'month');
biomass=ncread(filename,'biomass');

ny=length(lat);
ntim=length(month);
ndep=length(depth);

% what plankton groups to include
pnum=[1:31]; %phytoplankton (including mixotrophs)
nplank=length(pnum);
if (ntim==1 & ndep==1), 
 pbiomass=squeeze(biomass(pnum,:));
end % if
if (ntim>1 & ndep==1),
 pbiomass=squeeze(biomass(pnum,:,:));
end
if (ntim==1 & ndep>1),
 pbiomass=squeeze(biomass(pnum,:,:));
end
if (ntim>1 & ndep>1),
 pbiomass=squeeze(biomass(pnum,:,:,:));
end


% turn 3/4 dim matrixes into 2-dimensional matrix 
% maintain lat/depth/month info for each record
ivect=0;
for it=1:ntim
 for idep=1:ndep
  for il=1:ny,
   ivect=ivect+1;
   reclat(ivect)=lat(il);
   recdepth(ivect)=depth(idep);
   rectime(ivect)=month(it);
   for ip=1:nplank,
    clear tmp, 
    if (ntim==1 & ndep==1), tmp=squeeze(pbiomass(ip,il)); end
    if (ntim>1 & ndep==1), tmp=squeeze(pbiomass(ip,il,it)); end
    if (ntim==1 & ndep>1), tmp=squeeze(pbiomass(ip,il,idep)); end
    if (ntim>1 & ndep>1), tmp=squeeze(pbiomass(ip,il,idep,it)); end
    vect_biomass(ivect,ip)=tmp;
   end % for ip
  end
 end
end
vect_len=ivect;

% brey curtis similarity from each entry to each other
for il1=1:vect_len,
 tmp1=squeeze(vect_biomass(il1,:)); stmp1=sum(tmp1);
 for il2=1:vect_len,
   tmp2=squeeze(vect_biomass(il2,:)); stmp2=sum(tmp2);
   bctmp=0;
   for ip=1:nplank,
     bctmp=bctmp+min(tmp1(ip),tmp2(ip));
   end % for ip
   cb(il1,il2)=((2*bctmp)/(stmp1+stmp2));
 end % for il2
end % for il1

figure
pcolor(cb); shading flat, colorbar

% save into netcdf file to pass to infomap (python)
 outstr=['workshop_breycurtis_forinfomap.nc'];
 ncid=netcdf.create([outstr],'NC_WRITE');

 dimny=netcdf.defDim(ncid,'domainLat',ny);
 dimnz=netcdf.defDim(ncid,'domainDepth',ndep);
 dimnt=netcdf.defDim(ncid,'domainMonth',ntim);
 dimnp=netcdf.defDim(ncid,'nrecord',vect_len);

 lat_ID=netcdf.defVar(ncid,'domain_latitude','double',[dimny]);
 depth_ID=netcdf.defVar(ncid,'domain_depth','double',[dimnz]);
 time_ID=netcdf.defVar(ncid,'domain_month','double',[dimnt]);
 rec_ID=netcdf.defVar(ncid,'nrec','double',[dimnp]);

 cb_ID=netcdf.defVar(ncid,'CBmatrix','double',[dimnp,dimnp]);
 reclat_ID=netcdf.defVar(ncid,'matrixLat','double',[dimnp]);
 recdepth_ID=netcdf.defVar(ncid,'matrixDepth','double',[dimnp]);
 rectime_ID=netcdf.defVar(ncid,'matrixTime','double',[dimnp]);

 netcdf.endDef(ncid);

 netcdf.putVar(ncid,lat_ID,lat);
 netcdf.putVar(ncid,depth_ID,depth);
 netcdf.putVar(ncid,time_ID,month);
 netcdf.putVar(ncid,cb_ID,cb);
 netcdf.putVar(ncid,reclat_ID,reclat);
 netcdf.putVar(ncid,recdepth_ID,recdepth);
 netcdf.putVar(ncid,rectime_ID,rectime);

 netcdf.close(ncid);

