% session 4: more strings, parameters files
% and interpolation

clear all; clf

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 % interacting with parameter files
 
 
 dostart=1;
 valvec=floor(5:3:30);
 for ind=1:length(valvec);
    
     val=valvec(ind)
 
     if dostart
     st=['sed -e ''s/MYVALUE/' num2str(val) '/g'' datafile.i > fsc' num2str(ind) '.i']
 eval(['! ' st]);
 
 % starting the computation
 eval(['! fsc fsc' num2str(ind) '.i']);
 
     else
         % all the commands to get back the data to matlab
 
 
 end
 
 
 
 %%%%%%%%%%%%%%%%%%%
 %  example of manipulation folders with "!"
 break
 
 for ind=1:5
     repname=['myrep' num2str(ind)];
     com=['! mkdir ' repname ' ; cp fsc/fsc.i ' repname '/.']
     eval(com);
 end
 
     
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%ù
 break
 % simple example of building a command with strings

 st=[];
 for ind=1:30;
     st=[st 'disp(' num2str(ind) ');'];
 end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%
 break
 % example of interpolation: rotation
 
 a=imread('rayleigh.jpg');
 
 n=100
 for gre=1:20
 for phi=linspace(0,2*pi,50);;
 
 s=size(a);nl=s(1); nc=s(2);
 l=linspace(50,170,n);
 c=linspace(70,170,n);
 
 [C,L]=meshgrid(c,l);
 
 c0=105; l0=105;
 [TH,R]=cart2pol(C-c0,L-l0);
 [C,L]=pol2cart(TH+phi,R);
 C=C+c0; L=L+l0;
 
 z=interp2(1:nc,1:nl,double(a),C,L);
 
 subplot(1,2,2);
 imagesc(c,l,z);
 axis equal tight
 %pause(0.5)
 
 
 
 subplot(1,2,1);
 imagesc(a);hold on
 mesh(C,L,0*C);hold off
 
 colormap gray
 axis equal tight
 
 
 
 drawnow
 end
 
 end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 break
 % example of interpolation: pixelisation
 
 n=200;
 x=linspace(-2,2,n);
 y=linspace(-2,2,n);
 [X,Y]=meshgrid(x,y);
 f=cos(5*X).*exp(-X.^2-Y.^2);
 %surf(X,Y,f); 
 camlight; 
 
 for t=linspace(0,10*pi,500);
 clf; camlight
 
 th=linspace(0,15*2*pi,2000);
 r=(th/max(th))*2;
 x=r.*cos(th-t);
 y=r.*sin(th-t);
 %plot(x,y,'r-');
 
 [ff]=interp2(X,Y,f,x,y);
 
 hold off
 %plot3(x,y,ff,'r','linewidth',3);
 streamtube({[x',y',ff']},{x*0+0.1});
 shading interp; view(3); colormap gray
 axis equal tight
 colormap jet
 axis([-2,2 -2 2 -1,1])
 drawnow
 end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 break
 
 % 
 % 
 % break
 % 
 % 
 % 
 % % une fonction recursive pour créer des répertoires
 % 
 % 
 % %function st=createbranch(n,l,lmax,st)
 % % n the number of folders
 % % l the present level in the tree
 % % st the command string
 % 
 %     for ind=1:n
 % 
 %         % the name of the folder
 %         repname=['r_' num2str(l) '_' num2str(ind)];
 %         
 %         % create the folder
 %         st=[st 'mkdir ' repname ' ; '];
 %         
 %         % if not reached the max level
 %         % move to the new folder and create a branch
 %         % then move down
 %         if l<lmax
 %         st=[st ' cd ' repname ' ; '];
 %         st=createbranch(n,l+1,lmax,st);
 %         st=[st ' cd .. ; ']; 
 %         end
 %     end
 % 
 % 
 % 
 % 
 % 
 % 
 % 
 % 
 % break

 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % particules interpolation and advection

for gre=1:100
n=30;
x=linspace(0,2*pi,n);
y=linspace(0,pi,n/2);

[X,Y]=meshgrid(x,y);


% initial condition for particles
np=35000;
xp=rand(np,1)*2*pi;
yp=rand(np,1)*pi;

selb=xp<pi;

hold on

dt=0.1;

meth='linear';

nt=10000;
tvec=0:dt:nt*dt;
% my time loop
for ind=1:nt
t=tvec(ind);
u=sin(X-t).*cos(Y);
v=-cos(X-t).*sin(Y);
quiver(X,Y,u,v,1,'k');
hold on
    
up=interp2(X,Y,u,xp,yp,meth);
vp=interp2(X,Y,v,xp,yp,meth);



%quiver(xp,yp,up,vp,1,'r')

xp=mod(xp+dt*up,2*pi);
yp=yp+dt*vp;

plot(xp(selb),yp(selb),'b.','markersize',5);
%plot(xp(~selb),yp(~selb),'r.','markersize',5);
hold off
axis([0 2*pi 0 pi])
drawnow
%print(['frame' num2str(ind)],'-dpng');
end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
break
% example of interpolation in 2D

% un meilleur exemple de interpolation en ligne que avec Rayleigh
n=300;
x=linspace(-2,2,n);
y=linspace(-2,2,n);
[X,Y]=meshgrid(x,y);
f=cos(5*X).*exp(-X.^2-Y.^2);
%R=8*(X.^2+Y.^2); f=sin(R)./R;
%surf(X,Y,f); shading interp; camlight
colormap jet

camlight

for phi=linspace(0,2*pi,100);
th=linspace(0,30*pi,2000);
r=th/(30*pi)*2;
x0=0; y0=0;
x=r.*cos(th-phi)+x0;
y=r.*sin(th-phi)+y0;

z=interp2(X,Y,f,x,y);
hold on;
try; delete(h); end
%h=plot3(x,y,z,'r-','linewidth',2);
h=streamtube({[x',y',z']},{x*0+0.1});
shading interp; 
lighting gouraud
 set(gca,'cameraposition',[  6.4433  -33.1471    3.8649
]);
axis([-2,2,-2,2,-1,1])
box on
drawnow

end

