n = 15
m = 12;
r = linspace(0,1,n);

rrintr = 0;
ththintr = 0;
for k = 2:n-1
  c = (k-1)*m;
  th = [0:2*pi/c:2*pi*(1-1/c)];
  for j = 1:c
    rrintr(end+1) = r(k);
    ththintr(end+1) = th(j);
  endfor
endfor
#Remove one in order to get a nice non-prime number.
#rrintr = rrintr(2:end);
#ththintr = ththintr(2:end);

r=1;
c = (n-1)*m;
th = [0:2*pi/c:2*pi*(1-1/c)];
rrbdry = 1;
ththbdry = 0;
for j = 2:c
  rrbdry(end+1) = r;
  ththbdry(end+1) = th(j);
endfor
[xxintr,yyintr] = pol2cart(ththintr, rrintr);
[xxbdry,yybdry] = pol2cart(ththbdry, rrbdry);
intr = [xxintr',yyintr'];
bdry = [xxbdry',yybdry'];

##gotta avoid zero coordinates, so shift the boundary slightly
bdry = ([cos(pi/200), sin(pi/200); -sin(pi/200) cos(pi/200)]*bdry')';
nrml = [bdry, bdry];
save "circ_intr.matrix" intr;
save "circ_bdry.matrix" bdry;
save "circ_nrml.matrix" nrml;

xxbdry = bdry(:,1)';
yybdry = bdry(:,2)';

##radial
##u0 = (1-(xxintr.^2 + yyintr.^2)')/2; 

##nonradial
h0 = [];
 for i = 1:columns(xxintr)
   h0(i) = fr( (xxintr(i)-0.25)^2 + yyintr(i)^2, 0.025);
 endfor
 h0 = 4*h0(:)+0.1;

h0 = [h0; 0.1*ones(rows(bdry),1)];



h_init = [xxintr', yyintr'; xxbdry', yybdry'];
h_init = [h_init, h0];

save "h_init.map" h_init;


hold off;
plot(intr(:,1),intr(:,2),'xr');
 hold on;
plot(bdry(:,1),bdry(:,2),'xb');
size([intr;bdry])