kwr=1;
nc=24; 
np=nc+1; n6x=6*nc; nPg=n6x*nc;
addAngle=1;

outpName='grid_cs24';
inpName='tile.mitgrid';
rDir='./';

for n=1:6,
%-- read :
 p=strfind(inpName,'.mitgrid');
 if isempty(p),
   namF=sprintf([rDir,inpName,'.face%3.3i.bin'],n);
 else
   namF=sprintf([rDir,inpName(1:4),'%3.3i',inpName(5:end)],n);
 end
 fid=fopen(namF,'r','b');
 vv1=fread(fid,'real*8');
 fclose(fid);
 s=size(vv1,1); k1=s/np/np;
 fprintf(['read: ',namF,' : size: %i (%ix%ix%i)\n'],s,np,np,k1);
 vv1=reshape(vv1,[np np k1]);
%--- save each face:
 if n == 1,
   vvf=zeros(np,np,6,k1);
 end
 vvf(:,:,n,:)=vv1;
%-- fix missing corner values :
 if mod(n,2) == 1,  i=1; j=np; else i=np; j=1; end
%- xG k=  6 : Corner Max diff m1 = 1.42109e-14 (f3,f5= 1.42109e-14 -7.10543e-15)
 if n == 3 | n == 5,
   k=6;
   vvf(i,j,n,k)=vvf(i,j,1,k);
 end
%- yG k=  7 : Corner Max diff m2 = 7.10543e-15 (f4,f6= 0 7.10543e-15)
 if n == 6,
   k=7;
   vvf(i,j,n,k)=vvf(i,j,2,k);
 end
end

%- XC YC DXF DYF RA XG YG DXV DYU RAZ DXC DYC RAW RAS DXG DYG:
if addAngle == 1,
%-- Compute angle:
 var=vvf(1:nc,1:nc,:, 7); var=reshape(permute(var,[1 3 2]),[nPg 1]);
  yG =zeros(nPg+2,1); yG(1:nPg,1)=var;
  yG(nPg+1,1)=vvf(1,np,1,7);  yG(nPg+2,1)=vvf(np,1,2,7);
 var=vvf(1:nc,1:nc,:, 5); rAc=reshape(permute(var,[1 3 2]),[n6x nc]);
 var=vvf(1:nc,1:nc,:,15); dxG=reshape(permute(var,[1 3 2]),[n6x nc]);
 var=vvf(1:nc,1:nc,:,16); dyG=reshape(permute(var,[1 3 2]),[n6x nc]);

 [anCs,anSn] = cubeCalcAngle(yG,rAc,dxG,dyG);
 anCs=permute(reshape(anCs,[nc 6 nc]),[1 3 2]);
 anSn=permute(reshape(anSn,[nc 6 nc]),[1 3 2]);
end

if kwr == 1,
%-- write to file:
 for n=1:6,
  if addAngle == 1,
  %-- Add angle:
   k2=k1+2;
   vv2=zeros(np,np,k2); vv2(:,:,[1:k1])=vvf(:,:,n,:);
   vv2([1:nc],[1:nc],k1+1)=anCs(:,:,n);
   vv2([1:nc],[1:nc],k1+2)=anSn(:,:,n);
  else
  %-- do not add angle:
   vv2=vvf(:,:,n,:);
  end
  namW=sprintf([outpName,'.face%3.3i.bin'],n);
  fid=fopen(namW,'w','b');
  fwrite(fid,vv2,'real*8');
  fclose(fid);
  fprintf([' write to file: ',namW,' %i 2D.var(%ix%i)\n'],k2,np,np);
 end
end