function [salida]=cristales_puestos2definitivo(crystal);
%function [h1,h2,m,m_1,w,dt]=cristales(crystal)
%crystal=[1 1 1 1]. Se pone difusor. Si es 0 no hay difusor
%crystal=[2 2 2 2]. Se pone difusor optimizado
ZZ=60;
close all
if size(crystal,2)~=4
    disp('Error. Valor incorrecto del flag crystal. Debe ser una vector 1x4')
    return
end


maxtt=160000;%el numero de muestras tiene que ser superior a 81840
%para que dure mas de 0.6 segundos la simulación
pasovisionado=200;

Z=20;%Z= impendanca normalizada, tener en cuenta para que el RT de la sala sea
%el que se quiere, luego tneer en cuenta el número de muestras
maxx=3.7+.35*2;maxy=maxx; %DIMENSIONES DE LA SALA
sourcex=maxx-.75; %posición fuente x
sourcey=.75; % posición fuente y
dh=0.01; %PASO ESPACIAL
nx=round(maxx/dh);ny=round(maxy/dh);

rho=1.021;c=341;k=(c^2)*rho;dt=dh/c/2; % medio constantes
%initalization
p=zeros(nx,ny);
ux=zeros(nx+1,ny); %particle velocity
uy=zeros(nx,ny+1); %particle velocity

[m,m_1]=logsinesweep2(1/dt,(maxtt)/(1/dt));
t=(1:maxtt)*dt;
w=[m m*0];

figure(10);
subplot(2,1,1)
plot(t,w,t(pasovisionado:pasovisionado:end),w(pasovisionado:pasovisionado:end),'or');
subplot(2,1,2)
W=abs(fft(w)); n=length(W);
plot((0:n/2-1)/n/dt,10*log10(W(1:end/2).^2))
pause(2)
close




[x,y]=meshgrid((1:nx)*dh,(1:ny)*dh);x=x';y=y';
listap0=[];

a_red=.08;
radio=.035;

uxz=[];uyz=[]; %inicializacion de la lista de puntos de velocidad que son frontera
ley=(-22:22)+.2*(cumsum((0>sin((-22:22)/22*2*pi))-.5));
centrosx=[];centrosy=[];

%%
if crystal(1)>0; %difusor a la izquierda.
	for ii=1:4 
        for jj=1:45
            x0=-.5*a_red+a_red*ii;
            y0=maxy/2+a_red*ley(jj);
            if (x0<maxx-radio)&(x0>0+radio)&(y0<maxy-radio)&(y0>0+radio)
                centrosx=[centrosx x0];
                centrosy=[centrosy y0];
            end
        end
    end
end
%%


%%
if crystal(2)>0; %difusor a la derecha.
	for ii=1:4 
        for jj=1:45
            x0=maxx+.5*a_red-a_red*ii;
            y0=maxy/2-a_red*ley(jj);
            if (x0<maxx-radio)&(x0>0+radio)&(y0<maxy-radio)&(y0>0+radio)
                centrosx=[centrosx x0];
                centrosy=[centrosy y0];
            end
        end
    end
end
%%

%%
if crystal(3)>0; %difusor parte baja.
	for ii=1:4 
        for jj=1:45
            y0=-.5*a_red+a_red*ii;
            x0=maxx/2-a_red*ley(jj);
            if (x0<maxx-radio)&(x0>0+radio)&(y0<maxy-radio)&(y0>0+radio)
                centrosx=[centrosx x0];
                centrosy=[centrosy y0];
            end
        end
    end
end
%%
% 
% 
%%
if crystal(4)>0; %difusor parte alta
		for ii=1:4 
        for jj=1:45
            y0=maxy+.5*a_red-a_red*ii;
            x0=maxx/2-a_red*ley(jj);
            if (x0<maxx-radio)&(x0>0+radio)&(y0<maxy-radio)&(y0>0+radio)
                centrosx=[centrosx x0];
                centrosy=[centrosy y0];
            end
        end
    end
end
%%

    load radios
    radio=repmat(radios,1,4);
if crystal(1)==1
    radio=radio*0+1;
end


    for ii=1:length(centrosx)
       listap0=[listap0 find(sqrt((x-centrosx(ii)).^2+(y-centrosy(ii)).^2)<=radio(ii).*a_red/2*.9)'];
    end

%Busquedad de los puntos de contorno.
p=p;p(listap0)=1;
uxz=[];puxz=[];signoux=[];uyz=[];puyz=[];signouy=[];
%Busquedad de los puntos de contorno.
%Busquedad de los puntos de contorno de ux 
[i,j]=find(diff(p)==1);
puxz=[puxz (i+(j-1)*nx)'];
uxz=[uxz (i+1+(j-1)*(nx+1))'];
signoux=[signoux ones(length(i),1)'];
[i,j]=find(diff(p)==-1);
puxz=[puxz (1+i+(j-1)*nx)'];
uxz=[uxz (i+1+(j-1)*(nx+1))'];
signoux=[signoux -ones(length(i),1)'];

%Busquedad de los puntos de contorno de uy
[i,j]=find(diff(p')'==1);
puyz=[puyz (i+(j-1)*nx)'];
uyz=[uyz   (i+j*nx)'];
signouy=[signouy ones(length(i),1)'];
[i,j]=find(diff(p')'==-1);
puyz=[puyz (i+j*nx)'];
uyz=[uyz (i+j*nx)'];
signouy=[signouy -ones(length(i),1)'];
%%
p=p*0; %reseteo p para que no tenga la información de los cristales
%FIN CONDICIONES CONTORNO CRISTALES
%%


% %%%CABEZA
% x0=2;y0=1.8; %OJO ESTA ES LA POSICION DE LA CABEZA. MODIFICAR SI SE NECESITA MOVER AL ESPECTADOR
% radiocabeza=.10;
% listapcabeza=find(sqrt((x-x0).^2+(y-y0).^2)<=radiocabeza)';
% %Busquedad de los puntos de contorno.
% mip=p;mip(listapcabeza)=1;
% rightearx=x0/dh;
% righteary=find(mip(rightearx,:)>0,1,'first')-1;
% leftearx=x0/dh;
% lefteary=find(mip(rightearx,:)>0,1,'last')+1;
% 
% uxzcabeza=[];uyzcabeza=[];
% %Busquedad de los puntos de contorno de ux 
% [i,j]=find(diff(mip)==1);
% uxzcabeza=[uxzcabeza (i+1+(j-1)*(nx+1))'];
% [i,j]=find(diff(mip)==-1);
% uxzcabeza=[uxzcabeza (i+1+(j-1)*(nx+1))'];
% %Busquedad de los puntos de contorno de uy
% [i,j]=find(diff(mip')'==1);
% uyzcabeza=[uyzcabeza   (i+j*nx)'];
% [i,j]=find(diff(mip')'==-1);
% uyzcabeza=[uyzcabeza (i+j*nx)'];
% %Fin Busquedad de los puntos de contorno de cabeza.



for tt=1:maxtt;
%preassure calculation
   p=p-k*dt/dh*(diff(ux)+diff(uy')');
   p(round(sourcex/dh),round(sourcey/dh))=w(tt); %fuente posición
%velocity calculation
   ux(2:nx,:)=ux(2:nx,:)-dt/rho/dh*diff(p);
   uy(:,2:ny)=uy(:,2:ny)-dt/rho/dh*diff(p')';
%boundary conditions
   ux(1,:)=-p(1,:)/rho/c/Z;
   ux(end,:)=p(end,:)/rho/c/Z;
   uy(:,1)=-p(:,1)/rho/c/Z;
   uy(:,end)=p(:,end)/rho/c/Z;
%    ux(uxzcabeza)=0;   uy(uyzcabeza)=0;   CABEZA (DESCOMENTAR PARA AÑADIR)
   ux(uxz)=signoux.*p(puxz)/rho/c/ZZ; %IMPEDANCIA DEL CRISTAL=10
   uy(uyz)=signouy.*p(puyz)/rho/c/ZZ; %IMPEDANCIA DEL CRISTAL=10


%recording poner las coordenadas espaciales según el recinto
 
   h1(:,tt)=p(end/5,end/5); %La separación de los micrófonos debe tener en cuenta la posición y tamaño(30cm) de la cabeza
   h2(:,tt)=p(end/5,2*end/5);
   h3(:,tt)=p(end/5,3*end/5);
   h4(:,tt)=p(end/5,4*end/5);
   h5(:,tt)=p(2*end/5,end/5);
   h6(:,tt)=p(2*end/5,2*end/5);
   h7(:,tt)=p(2*end/5,3*end/5);
   h8(:,tt)=p(2*end/5,4*end/5);
   h9(:,tt)=p(3*end/5,end/5);
   h10(:,tt)=p(3*end/5,2*end/5);
   h11(:,tt)=p(3*end/5,3*end/5);
   h12(:,tt)=p(3*end/5,4*end/5);
   h13(:,tt)=p(4*end/5,end/5);
   h14(:,tt)=p(4*end/5,2*end/5);
   h15(:,tt)=p(4*end/5,3*end/5);
   h16(:,tt)=p(4*end/5,4*end/5);



if tt/pasovisionado==round(tt/pasovisionado)
    warning off;
    pp=p;
    pp=max(10*log10(pp.^2),-80);
    pp(listap0)=-120;
    pcolor((1:nx)*dh,(1:ny)*dh,pp');
    set(gca,'Clim',[-120 0]);
    shading flat,axis equal; colorbar;
    title(['time=' num2str(round(tt*dt*1000)) ' miliseconds' '   (listo ' num2str(tt) ' de '  num2str(maxtt)  ')']);drawnow
    warning on
end


end

% figure;
% subplot(2,1,1);plot((1:length(h1))*dt,h1); title('Respuesta en el tiempo canal R')
% subplot(2,1,2);plot((1:length(h2))*dt,h2); title('Respuesta en el tiempo canal L')


mic1=conv(h1,m_1);
mic2=conv(h2,m_1);
mic3=conv(h3,m_1);
mic4=conv(h4,m_1);
mic5=conv(h5,m_1);
mic6=conv(h6,m_1);
mic7=conv(h7,m_1);
mic8=conv(h8,m_1);
mic9=conv(h9,m_1);
mic10=conv(h10,m_1);
mic11=conv(h11,m_1);
mic12=conv(h12,m_1);
mic13=conv(h13,m_1);
mic14=conv(h14,m_1);
mic15=conv(h15,m_1);
mic16=conv(h16,m_1);

normalizacion=max([mic1 mic2 mic3 mic4 mic5 mic6 mic7 mic8 mic9 mic10 mic11 mic12 mic13 mic14 mic15 mic16]);
mic1=mic1(1:end)/normalizacion*.9;
mic2=mic2(1:end)/normalizacion*.9;
mic3=mic3(1:end)/normalizacion*.9;
mic4=mic4(1:end)/normalizacion*.9;
mic5=mic5(1:end)/normalizacion*.9;
mic6=mic6(1:end)/normalizacion*.9;
mic7=mic7(1:end)/normalizacion*.9;
mic8=mic8(1:end)/normalizacion*.9;
mic9=mic9(1:end)/normalizacion*.9;
mic10=mic10(1:end)/normalizacion*.9;
mic11=mic11(1:end)/normalizacion*.9;
mic12=mic12(1:end)/normalizacion*.9;
mic13=mic13(1:end)/normalizacion*.9;
mic14=mic14(1:end)/normalizacion*.9;
mic15=mic15(1:end)/normalizacion*.9;
mic16=mic16(1:end)/normalizacion*.9;

mic1=resample(mic1,44100,1/dt);
mic2=resample(mic2,44100,1/dt);
mic3=resample(mic3,44100,1/dt);
mic4=resample(mic4,44100,1/dt);
mic5=resample(mic5,44100,1/dt);
mic6=resample(mic6,44100,1/dt);
mic7=resample(mic7,44100,1/dt);
mic8=resample(mic8,44100,1/dt);
mic9=resample(mic9,44100,1/dt);
mic10=resample(mic10,44100,1/dt);
mic11=resample(mic11,44100,1/dt);
mic12=resample(mic12,44100,1/dt);
mic13=resample(mic13,44100,1/dt);
mic14=resample(mic14,44100,1/dt);
mic15=resample(mic15,44100,1/dt);
mic16=resample(mic16,44100,1/dt);

salida=[mic1', mic2', mic3', mic4', mic5', mic6', mic7', mic8', mic9', mic10', mic11', mic12', mic13', mic14', mic15', mic16'];


% canalL=conv(h1,m_1);
% canalR=conv(h2,m_1);
% normalizacion=max([canalL canalR]);
% canalL=canalL(1:end)/normalizacion*.9;
% canalR=canalR(1:end)/normalizacion*.9;
% canalL=resample(canalL,44100,1/dt);
% canalR=resample(canalR,44100,1/dt);
% salida=[canalL',canalR'];
%wavwrite(salida,44100,'resultado.wav')



return

