% Field map with interpolation, rods and jacket position. 
% For an "easy use", only change the line giving the titre, just below.
% This script works with Matlab 6.1 and 7.6
%=====================================================================
clear all; close all; format compact;
titre = 'fig11_JOSAA_2004'; % Put here the string used in the first line of the data input file
rodsfile = [titre,'.xy.dat'];
rods = dlmread(rodsfile,'',1,0);     % use ' ' for Matlab 6.1, and '' for Matlab 7.6 
fieldfile = ['cartechamp.',titre,'.dat'];
jacketfile = ['prof_norm.',titre];
field = dlmread(fieldfile,'',1,0);   % use ' ' for Matlab 6.1, and '' for Matlab 7.6 
jacket = dlmread(jacketfile,'',1,0); % use ' ' for Matlab 6.1, and '' for Matlab 7.6 

% Rods position :
xc = rods(:,2);
yc = rods(:,3);
radius = rods(:,9); dia = 2*radius;
nrods=length(xc);

% Jacket profile :
xd = jacket(:,1); yd = jacket(:,2);

% Field
x = unique(field(:,1)); xmin = min(x); xmax = max(x); nx=length(x); dx = (xmax-xmin)/(nx-1);
y = unique(field(:,2)); ymin = min(y); ymax = max(y); ny=length(y); dy = (ymax-ymin)/(ny-1);
modu=field(:,3); % 3 plots the modulus. 4 plots the phase
z=reshape(modu,ny,[]);
% We interpolate on a grid 2 times thiner
[xi,yi] = meshgrid(xmin:dx/2:xmax,ymin:dy/2:ymax);
zi = interp2(x,y,z,xi,yi);

% Choose between the two following graphs:

% 1) Normal graph :
pcolor(xi,yi,zi) 
shading interp 
caxis([0 1]); %Optional: choose limits of the colormap
colorbar;

% 2) Contour map
%zmin = 0; zmax = 1; deltaz = 0.1; % used for isolines and colorbar
%nbcolors = (zmax-zmin)/deltaz;
%colormap(jet(nbcolors)); 
%v = [zmin:deltaz:zmax]  
%[C,h] = contourf(xi,yi,zi,v); % contour map with isolines
%clabel(C,h); Optional: add the isolines labels
%caxis([zmin zmax]); colorbar;

axis square; 
title(titre,'Interpreter','none');

% Rods
for it=1:nrods
    rectangle('Position',[xc(it)-radius(it),yc(it)-radius(it),dia(it),dia(it)],'Curvature',[1,1],'EdgeColor','k','LineWidth',0.5)
end

% Jacket
hold on;
plot(xd,yd,'k','LineWidth',0.5);

%----------------------------------------------------
% Place des tickmarks sur les 4 cotés, et les met en avant-plan:
set(gca,'Box','On','Layer','top');

%----------------------------------------------------
% Exporter en fichier graphique:
% On peut le faire depuis la figure avec File / Export
% ou bien en ligne de commande:
print('-depsc','-r300',[titre,'-epsc_300']) % eps 300 dpi
print('-dpng','-r300',[titre,'-png_300']) % png 300 dpi

