% circ2.m
% Creates a circular mesh 
% Now with meshcd  Drag coefficient included

dx = input('Length of side of one element (750) >> ')
dy = dx;
lx = input('Width of the inner fine grid (23000) >>')
ly = lx;

[x,y,meshele,meshnod]=trigrid(dx,dy,lx,ly,30);
[meshele,meshnod]=cleanup(meshele,meshnod);


% Find the boundary info
 [bndele] = findbndy(meshele,meshnod);
figure(1)
plotmesh
figure(2)
[ax] = ploteles(meshele,meshnod)

Xc = mean(meshnod(:,2));
Yc = mean(meshnod(:,3));
meshnod(:,2) = meshnod(:,2)-Xc;
meshnod(:,3) = meshnod(:,3)-Yc;
Radius = ( max(meshnod(:,2)) - min(meshnod(:,2)) ) /2
figure(3)
[meshnod]= pullnods(meshele,meshnod,bndele,Radius,0,0);

[meshnod]=sprngnod(meshele,meshnod,bndele,30);


disp([' Radius = ',num2str(Radius),...
     '   number boundary points = ',int2str(length(bndele))])
numcircles = input('Number of circles outside submesh >>');
numradii = input('Number of submesh radii for outermost circle >>');
numbndpts = input('Number of boundary points for outermost circle >>');

nbnd = (length(bndele(:,1))-1)
fact = exp(log(numbndpts/nbnd)/(numcircles))
npts = fix(nbnd*fact.^([1:numcircles]))
newnodes = ones(sum(npts),3);
the(1)=0.;
for ic = 1:numcircles
r = Radius*(1+ic*numradii/numcircles);
r = Radius*(1+numradii*sum(npts((numcircles-ic+1):numcircles))/sum(npts));
the = 2*pi*[1:npts(ic)]/npts(ic)+the(1)/2;
lnode=[(sum(npts(1:ic))+1-npts(ic)) : sum(npts(1:ic))]';
newnodes(lnode,2) = r*cos(the');
newnodes(lnode,3) = r*sin(the');
end
plot(newnodes(:,2),newnodes(:,3),'ro',meshnod(:,2),meshnod(:,3),'bo')
meshnod = [meshnod; newnodes];

%dirnow = cd;
%tic
%cd /home5/tom/gridmaker
%[ meshele, meshnod, bndele]=gridmake(meshnod);
%eval(['cd ',dirnow]); toc


% USE  SAGA  TO DO TRIANGULATION

[meshele] = makeeles(meshnod);
[bndele]= findbndy(meshele,meshnod);




nrelaxes = input('How many times to relax the grid? >>');
figure(2)
[meshnod]=sprngnod(meshele,meshnod,bndele,nrelaxes);

% Initialize at constant depth
depth = 10;
bathy = [[1:length(meshnod)]' ones(length(meshnod),1)* depth];


% Create two longitudinal exponential shaped humps
% Use with dx = 500   L = 80*500 = 40000
%humplength = 12000;
%humpwidth = 2000;
%humpspace =6000;
%humpheight = 10;
%depth = 20;
%[bathy]=bumpexp(meshnod,depth,humpheight,humplength,humpspace,humpwidth);


% COSINE BUMPS
depth = 20;
humpheight = 10;
humplength = 12000;
humpspace =6000;
humpwidth = 3000;
depth = input(' What is depth of the bottom (20)>>');
humpheight = input(' What is height of the bumps above the bottom (10)>>');
cdd = input(' What is far field drag coefficient (0.005) >>');
cddbump = input(' What drag coefficient of the bumps (0.0025, )>>');

[bathy]=bumpcos(meshnod,depth,humpheight,humplength,humpspace,humpwidth);
[meshcd]=bumpcos(meshnod,cdd,cddbump,humplength,humpspace,humpwidth);

figure(3)
clf
subplot(1,2,1)
plot3(meshnod(:,2),meshnod(:,3),-bathy(:,2),'o'); title('depth');axis('equal')
grid; xlabel(cd);ylabel(date)
subplot(1,2,2)
plot3(meshnod(:,2),meshnod(:,3),meshcd(:,2),'o'); title('Cd');axis('equal')
grid

figure(4)
 [bndele] = findbndy(meshele,meshnod);
% amplitude and phase of sea surface for M2 forcing
bcs = ampphase(bndele,meshnod,bathy);


disp(' cleanup and bandwidth tighten')

% remove unwanted triangles and nodes
[meshele,meshnod,bathy,meshcd,bndele]=...
     reelems(meshele,meshnod,bathy,meshcd,bndele);

figure(1)
% Tighten up the bandwidth of the mesh
[meshele,meshnod,bathy,meshcd,bndele] = ...
        reorder(meshele,meshnod,bathy,meshcd,bndele);

 [bndele] = findbndy(meshele,meshnod);
% amplitude and phase of sea surface for M2 forcing
bcs = ampphase(bndele,meshnod,bathy);

% call the file writing routine
filename = input('What is root name of the mesh files >> ','s')
fileout

figure(1)
plot(meshnod(:,2),meshnod(:,3),'o'); axis('equal')
L = axis;
[xc, yc]=contmesh(L,meshele,meshnod,bathy(:,2),[0:2:18 19]);
axis('equal')
colormap(pink)
title(filename)
figure(2)
plot(meshnod(:,2),meshnod(:,3),'y.',xc',yc','r'); axis('equal')
title(filename)
figure(4)
[xc, yc]=contmesh(L,meshele,meshnod,meshcd(:,2),[0.001:0.001:0.01]);
axis('equal')
colormap(hot)
title(['Cd variability   ',cd,'/',filename,'.cd'])
