%Algoritmo para calcular la curva media extrinseca de un conjunto de curvas dado.

function curva2 %Funcion principal

%En esta funcion se piden los datos iniciales n y N, se crea el vector
%de angulos equiespaciados, t, con distancia h y se llama, secuencialmente,
%a las funciones que calculan los datos que necesitamos para hallar
%la curva media. Estas funciones estan descritas despues de esta funcion,
%en el mismo orden en que son llamadas.
%Ademas, abrimos un fichero de texto para guardar los datos que vamos
%obteniendo y creamos una grafica donde se representan todas las curvas;
%siendo la curva media, la negra.

n=input('Introduzca el número de curvas: ');
N=input('Introduzca el número de puntos que se toman en las curvas, teniendo en cuenta que se debe cerrar el recorrido: ');
h=input('Introduzca el valor del tamaño de paso h: ');
%distancia entre los angulos del intervalo [0,2pi]

%h=2*pi/(N-1); 
%En el ejercicio analizado hay varios puntos de las curvas con el mismo angulo 
%por lo que no se puede tomar el tamaño de paso con esta formula, porque el 
%numero de puntos final es 33 y se han tomado con tamaño de paso 2pi/28. 

filename=input('Introduzca el nombre del fichero excel que contiene los datos:  ','s');

t=linspace(0,2*pi,N); 
%Vector de angulos tomados con distancia h. Es importante entender este
%vector como uno de posiciones, es decir, el t(1)=0 y el t(N)=2pi.

fich=fopen('datos','w+');

c=curvas(n,N,t',filename);   
%Funcion para crear el conjunto de curvas. Desde la linea 256 hasta la
%312.

%- - - - - - - - - - - - - - Abrimos el fichero y comenzamos a guardar datos - - - - - - - - - - - - - %
fprintf(fich, ' - - - - - - RESULTADOS - - - - - - \n');
fprintf(fich,'\nDatos iniciales: \n');
fprintf(fich,'\nNúmero de curvas, número de puntos y tamaño de paso: %d, %d y %f.\n\n', n, N, h);
fprintf(fich, '- - - - - - - - - - - - - - - - - - - \n\n');
fprintf(fich, 'Curvas que creamos con esos puntos para este ejemplo.\n\n');

for i=1:n
    
    fprintf(fich, 'Curva %d\n\n', i);
    
    for j=1:N
        
        fprintf(fich, '| %f, %f |\n', c(j,1,i),c(j,2,i));
        
    end
    
    fprintf(fich, '\n');
 
end
fprintf(fich, '- - - - - - - - - - - - - - - - - - - \n\n');

fprintf(fich, 'Datos obtenidos por el algoritmo:\n\n');

%- - - - - - - - - - - - - - Continuamos llamando a las funciones - - - - - - - - - - - - - %

[d,theta]=deriv(fich,c,N,n,h);  
%Funcion para calcular la norma de derivada de las curvas y el angulo theta 
%que usaremos para obtener el par (a,b) correspondiente a cada curva. 
%De la linea 316 a la linea 470.
                                
L=longitud(h,n,d)';  
%Funcion para calcular la longitud de cada curva. De la linea 474 a la
%linea 494.

cent=centroide(c,L,d,N,n,h);  
%Funcion para calcular el centroide de cada curva. De la linea 498 a la
%linea 533.

[a,b]=vectorayb(d,L,theta,N,n);   
%Funcion para calcular el par (a,b) correspondiente a cada curva. 
%De la linea 537 a la linea 565.

[cm, Lm, am, bm]=cmedia(cent,L,a,b,n,N);   
%Funcion para calcular las medias directas del centroide y del par (a,b).
%De la linea 569 a la linea 593.

[a0m,b0m]=nuevosayb(am,bm,N,h);   
%Funcion para calcular la rotacion del par (am,bm). De la linea 597 a la
%linea 682.

curvaphi2=curvaporphi2(a0m,b0m,N,h);   
%Funcion para calcular la curva imagen por la aplicacion phi2. 
%De la linea 686 a la linea 719.

d2=der2(curvaphi2,N,h); 
%Funcion para calcular la norma de derivada de la curva imagen por phi2. 
%De la linea 723 a la linea 749.

L2=long(d2,h); 
%Funcion para calcular la longitud de la curva imagen por phi2. 
%De la linea 753 a la linea 761.

cent2=centroi(curvaphi2,L2,d2,N,h);  
%Funcion para calcular el centroide de la curva imagen por phi2. 
%De la linea 765 a la linea 788.

curvamedia=curvanueva(cm,Lm,curvaphi2,cent2,N);  
%Funcion para calcular la curva media. De la linea 792 a la linea 815.

d3=der3(curvamedia,N,h); 
%Funcion para calcular la derivada de la curva media. De la linea 832 a la
%linea 858.

L3=long(d3,h); 
%Funcion para calcular la longitud de la curva imagen. Comprobacion para este ejemplo. 
%De la linea 753 a la linea 761.

cent3=centroi(curvamedia', L3,d3, N, h);
%Function para calcular el centroide de la curva media. De la linea 765 a
%la linea 788.
r=zeros(n+1,1); %Radio de la curva media. Comprobacion de este ejemplo.

for i=1:n
    r(i)=L(i)/(2*pi);
end
r(n+1)=L3/(2*pi);

%- - - - - - - - - - - - - - Continuamos escribiendo en el fichero - - - - - - - - - - - - - %

fprintf(fich,'Norma de la derivada de cada curva.\n\n'); 

for i=1:n
    
    fprintf(fich,'Norma de la derivada de la curva %d\n\n',i);
    
    fprintf(fich,'(%f, ', d(1,i));
    for k=2:N-1
        
        fprintf(fich,'%f, ', d(k,i));
        
    end
    fprintf(fich,'%f)\n\n', d(N,i));
    
end
fprintf(fich, '- - - - - - - - - - - - - - - - - - - \n\n');

fprintf(fich, 'Ángulos theta que necesitamos para evaluar las curvas a y b correspondientes a cada curva.\n\n');

for i=1:n
    
    fprintf(fich, 'Ángulos theta correspondientes a la curva %d\n\n', i);
    fprintf(fich,'(%f, ', theta(1,i));
    for j=2:N-1
        
        fprintf(fich, '%f, ', theta(j,i));
        
    end
    fprintf(fich,'%f)\n\n', theta(N,i));
    
end

fprintf(fich, '- - - - - - - - - - - - - - - - - - - \n\n'); 

for i=1:n
    
    fprintf(fich, 'Longitud, radio y centroide de la curva %d: %f, %f y (%f,%f). \n\n', i, L(i),r(i),cent(1,i), cent(2,i));
    
end

fprintf(fich, '- - - - - - - - - - - - - - - - - - - \n\n'); 

for i=1:n
    
    fprintf(fich,'Par (a,b) correspondiente a la curva %d\n\n', i);
    
    for j=1:N
        
        fprintf( fich, '| %f, %f |\n', a(j,i), b(j,i));
        
    end
    fprintf(fich, '\n');
 
end

fprintf(fich, '- - - - - - - - - - - - - - - - - - - \n\n');

fprintf(fich, 'Longitud media y Centroide medio de las tres curvas (a,b): %f y (%f, %f)\n\n', Lm, cm(1,1), cm(2,1));

fprintf(fich, '- - - - - - - - - - - - - - - - - - - \n\n');

fprintf(fich, 'Par de vectores (am,bm)\n\n');

for i=1:N
    
    fprintf(fich, '| %f, %f |\n', am(i,1), bm(i,1));
    
end

fprintf(fich, '\n- - - - - - - - - - - - - - - - - - - \n\n');

fprintf(fich, 'Par de vectores (a0m,b0m)\n\n');

for i=1:N
    
    fprintf(fich, '| %f, %f |\n', a0m(i,1), b0m(i,1));
    
end

fprintf(fich, '\n- - - - - - - - - - - - - - - - - - - \n\n');

fprintf(fich, 'Curva phi2(a0m, b0m)\n\n');

for i=1:N
    
    fprintf(fich, '| %f, %f |\n', curvaphi2(1,i), curvaphi2(2,i));
    
end
fprintf(fich, '\n- - - - - - - - - - - - - - - - - - - \n\n');

fprintf(fich, 'Norma de la derivada de la curva phi2(a0m,b0m)\n\n');

fprintf(fich, '(%f, ', d2(1,1));
for i=2:N-1
    
    fprintf(fich, '%f, ', d2(i,1));
    
end
fprintf(fich, '%f)\n\n',d2(N,1));

fprintf(fich, '- - - - - - - - - - - - - - - - - - - \n\n');

fprintf(fich, 'Longitud y centroide de la curva phi2(a0m,b0m): %f y (%f, %f)\n\n',L2, cent2(1,1), cent2(2,1));

fprintf(fich, '- - - - - - - - - - - - - - - - - - - \n\n');

fprintf(fich, 'Curva media extrínseca buscada: \n\n');

for i=1:N
    
    fprintf(fich, '| %f, %f |\n', curvamedia(i,1), curvamedia(i,2));
    
end
fprintf(fich, '\n- - - - - - - - - - - - - - - - - - - \n\n');

fprintf(fich, 'Longitud, radio y centroide de la curva media extrínseca: %f, %f y (%f, %f).\n', L3, r(n+1), cent3(1,1), cent3(2,1));


%- - - - - - - - - - - - - - Creamos las gráficas - - - - - - - - - - - - - %

figure(1) 
%Creamos la grafica para representar las curvas junto con la curva media en color negro.

for i=1:n
    h=convhull(c(:,1,i),c(:,2,i)); %Para poder pintar una de las curva, debemos creas
                                   %los segmentos que unen los puntos de la misma.  
                                   %Para ello, usamos la envolvente convexa.
    colores=['y','r','b'];                               
    plot(c(h,1,i),c(h,2,i), colores(i)); % Con el comando plot se dibuja la grafica.
    hold on
end

c2nueva=curvamedia;
g=convhull(c2nueva(:,1),c2nueva(:,2));
plot(c2nueva(g,1),c2nueva(g,2),'k');  
hold off
title('Conjunto de curvas y curva media');
print('-djpeg', 'curva.jpg');

figure(2) 
for i=1:n
    h=convhull(c(:,1,i),c(:,2,i)); 
    colores=['y','r','b'];
                                   
    plot(c(h,1,i),c(h,2,i), colores(i));
    hold on
end
title('Conjunto de curvas');
print('-djpeg', 'curvasmuestra.jpg');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function c=curvas(n,N,t,filename)   % Funcion para crear el conjunto de curvas. 

%Datos que necesitamos: - numero de curvas que crearemos (n)
%                       - vector de angulos donde tomaremos los puntos de la
%                       curva (t)
%Salida: una matriz de matrices de dimensiones (Nx2)xn. Por ejemplo, el
%dato c(i, 1, k) es la componente x en la fila i, es decir, evaluada en t(i) de la curva k.

p=input('Marque un 1 si quiere ver el ejemplo de las circunferencias concéntricas o 2 si quiere introducir los puntos manualmente: %d\n'); 

if p==2
    
  v=zeros(N,2*n);
  %Creamos una matriz Nx2n para cargar los datos en filas, es decir, todas
  %las x la primera curva en la primera fila, todas las y en la segunda;
  %todas las x de la segunda curva en la tercera fila, las y en la cuarta;
  %y asi sucesivamente. Posteriormente, reorganizaremos los datos para
  %poder aplicar el algoritmo.

  for i=1:n
     
    %Para leer correctamente los puntos del excel, se deben colocar los
    %datos de estas en distintas hojas y separando las x de las y en
    %distintas columnas. 
    
      for j=1:N
        
          v(:,2*i-1)=xlsread(filename,i,'A1:A33');
          v(:,2*i)=xlsread(filename,i,'B1:B33');
         
          %Asegurese de indicar el rango correcto de celdas del excel.
      
      end
  end
  c=zeros(N,2,n);

  for i=1:n
      for j=1:N
        
          c(j,1,i)=v(j,2*i-1);
          c(j,2,i)=v(j,2*i);
        
      end
  end    
end

if p==1
    
    v=[cos(t),sin(t)]; %Creamos la parametrizacion

    c=cat(n,v);         %Primera curva. Se almacenara como la primera matriz (N,2)
    c=cat(n,c,2*v);     %Segunda curva. Se almacenara como la segunda matriz (N,2)
    c=cat(n,c,6*v);
    
    fprintf(fich,'\nVector de ángulos para evaluar las curvas y obtener los puntos ');
    fprintf(fich,'(%f, ', t(1));
    for i=2:N-1
        fprintf(fich,'%f, ', t(i));
    end
    fprintf(fich,'%f)\n\n', t(N));
    fprintf(fich, '- - - - - - - - - - - - - - - - - - - \n\n');

    
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [d,theta]=deriv(fich,c,N,n,h)  % Funcion para calcular la norma de la derivada de cada curva.

%Datos de entrada: - matriz de matrices de las curvas.
%                  - numero de angulos que se toman del intervalo [0,2pi]
%                  - numero de curvas
%                  - tamanio de paso h
%Salida: - matriz de dimensiones Nxn donde se almacena el valor de la norma
%          de la derivada de cada curva
%        - vector con los angulos theta necesarios para crear los pares
%          (a,b)

d=zeros(N,n);  
%Iniciamos la matriz de dimensiones Nxn, donde almacenaremos el valor de 
%la norma de la derivada de cada curva

v=zeros(N,2,n); 
%Iniciamos la matriz de matrices de dimensiones (Nx2)xn, donde guadaremos
%el valor de la derivada de cada curva. Esta matriz funciona igual que la
%matriz del conjunto de curvas, esto es, el dato v(i,2,k) es la y' en la
%fila i de la curva k.

%Como las curvas estan discretizadas, no podemos evaluar la derivada 
%forma continua, sino que debemos usar un metodo numerico. En este caso, 
%usaremos el siguiente metodo de orden cinco:

%(-c(i+2,:,k) + 8*c(i+1,:,k) - 8*c(i-1,:,k) + c(i-2,:,k))/12*h

%Destacar que al poner los dos puntos, la formula se aplica a la vez para
%las componentes x y las y.

%Sin embargo, al aplicar este metodo debemos tener cuidado con los
%extremos; pues al evaluar v(N,:,k), necesitamos el v(N+2, :, k). Como 
%las curvas con cerradas y c(1,:,k)=c(N,:,k), ese dato se corresponde con
%v(3,:,k), pero matlab nos daria un error porque entiende que ese dato no 
%existe. Por tanto, v(1,:,k), v(2,:,k),v(N-1,:,k) y v(N,:,k) hay que
%calcularlos de otra  forma.


for k=1:n
    for i=3:N-2
        
        %Con este for calculamos los que no dan problemas.
  
        v(i,:,k)= (-c(i+2,:,k) + 8*c(i+1,:,k) -8*c(i-1,:,k) +c(i-2,:,k) )/(12*h); 
               
        d(i,k)=sqrt(v(i,1,k)^2 + v(i,2,k)^2);  
        
    end
    
    %Y ahora calculamos los extremos.
    
    v(1,:,k)=( -c(3,:,k) + 8*c(2,:,k) -8*c(N-1,:,k) + c(N-2,:,k) )/(12*h);
    v(2,:,k)=( -c(4,:,k) +8*c(3,:,k) -8*c(1,:,k) +c(N-1,:,k) )/(12*h);
    
    d(1,k)= sqrt( v(1,1,k)^2 + v(1,2,k)^2);  
    d(2,k)= sqrt( v(2,1,k)^2 + v(2,2,k)^2);
    
    v(N-1,:,k)=( -c(2,:,k) +8*c(N,:,k) -8*c(N-2,:,k) + c(N-3,:,k) )/(12*h);
    v(N,:,k)=( -c(3,:,k) +8*c(2,:,k) -8*c(N-1,:,k) + c(N-2,:,k) )/(12*h);
    
    d(N-1,k)= sqrt( v(N-1,1,k)^2 + v(N-1,2,k)^2);
    d(N,k)= sqrt( v(N,1,k)^2 + v(N,2,k)^2);
    
end 

for i=1:n
    
    fprintf(fich, 'Derivada de la curva %d\n\n', i);
    
    for j=1:N
        
        fprintf(fich, '| %f, %f |\n', v(j,1,i),v(j,2,i));
        
    end
    fprintf(fich, '\n');
 
end
fprintf(fich, '- - - - - - - - - - - - - - - - - - - - -\n\n');

%A continuacion, calculamos los angulos theta que necesitamos para obtener
%los pares (a,b). Usaremos la funcion arcotangente(y'/x'), teniendo en
%cuenta que estamos trabajando en el intervalo [0,2pi] y no en
%[-pi/2,pi/2].
%Se debe tener en cuenta tambien que el ultimo valor de la curva se 
%corresponde con el primero, por ser una curva cerrada.

theta=zeros(N,n);   

for j=1:n
    
    for i=1:N
            
        if v(i,2,j) < 0  %Si y'<0
            
            if v(i,1,j) > 0  %Si x'>0
                
                theta(i,j)= 2*pi + atan(v(i,2,j)/v(i,1,j));
            end 
            
            if v(i,1,j) < 0 %Si x'<0
                
                theta(i,j)=pi + atan(v(i,2,j)/v(i,1,j));
            end
            
            if v(i,1,j)==0  %Si x'=0
                
                theta(i,j)=3*pi/2;
                
            end
        end
              
        if v(i,2,j) > 0 %Si y'>0
            
            if v(i,1,j) > 0  %Si x'>0
                
                theta(i,j)= atan(v(i,2,j)/v(i,1,j));
            end 
            
            if v(i,1,j) < 0 %Si x'<0
                
                theta(i,j)=pi + atan(v(i,2,j)/v(i,1,j));
            end
            
            if v(i,1,j)==0 %Si x'=0
                
                theta(i,j)=pi/2;
                
            end
            
        end
            
        if  v(i,2,j)==0 %Si y'=0
            
            if v(i,1,j) > 0 %Si x'>0
            
                theta(i,j)=atan(v(i,2,j)/v(i,1,j));
                %En este caso, como la curva es cerrada puede ser 0 o 2pi
                
            end
            
            if v(i,1,j) < 0 %Si x'<0
                
                theta(i,j)=pi + atan(v(i,2,j)/v(i,1,j));
                %=pi
                
            end
            
        end
        
    end 
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function L=longitud(h,n,d)  % Funcion para calcular la longitud de todas las curvas. 

%Datos de entrada: - numero de angulos que se toman del intervalo [0,2pi]
%                  - numero de curvas 
%                  - matriz de matrices de las curvas
%Salida: - vector de dimension n que contiene las longitudes de cada curva 

%Como no podemos calcular las integrales de manera continua, necesitamos
%introducir un metodo numerico para resolverlas. En este caso, hemos
%escogido el metodo de la Trapezoidal compuesto, computado de la linea a la 
%linea 

L=zeros(n,1);

for i=1:n
    
    L(i)=integral(h,d(:,i));
       
end    

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function cent=centroide(c,L,d,N,n,h)  %Funcion para calcular el centroide de cada curva.

%Datos de entrada: - matriz de matrices de curvas
%                  - vector de las longitudes de las curvas
%                  - matriz con la normas de las derivadas de las curvas
%                  - numero de angulos que se toman del intervalo [0,2pi]
%                  - numero de curvas
%                  - tamaño de paso h
%Salida: matriz de dimensiones 2xn que almacena el centroide de cada curva

%Teniendo en cuenta que el centroide es la integral de ||c'(t)||*c(t), la
%funcion que evalua la Regla del Trapecio compuesta sera un vector que almacene 
%el producto c(i,:,k)*d(i,k). 

cent=zeros(2,n);  

sum1=zeros(N,n);    
sum2=zeros(N,n);

%En estas matrices se almacenan las multiplicaciones c(i,:,k)*d(i,k) por
%columnas, correspondiendo cada columna con cada curva.

for k=1:n    
    for i=1:N    
        
        sum1(i,k)= d(i,k)*c(i,1,k); 
        sum2(i,k)= d(i,k)*c(i,2,k);   
         
    end
    
    cent(1,k)=integral(h,sum1(:,k))/L(k);    %Primera componente del centroide.
    cent(2,k)=integral(h,sum2(:,k))/L(k);    %Segunda componente del centroide.
    
end
end  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [a,b]=vectorayb(d,L,theta,N,n) %Funcion para calcular los pares de (a,b) correspondientes de cada curva

%Datos de entrada: - matriz con la norma de la derivada de cada curva
%                  - vector con la longitud de cada curva
%                  - vector de angulos theta que calculamos en la funcion
%                    [d,theta]
%                  - numero de angulos que tomamos del intervalo [0,2pi]
%                  - numero de curvas
%Salida: par (a,b)

%Destacar que, aunque teoricamente se trata del par (a,b), para matlab a y b
%son elementos diferenciados.

a=zeros(N,n);
b=zeros(N,n);

%En estas matrices se almacenan los pares (a,b) correspondientes a cada
%curva por columnas, es decir, el par (a(:,k),b(:,k)) es el par (a,b)
%completo relacionado con la curva k.
                              
for k=1:n
    for i=1:N
        g=2*d(i,k)/L(k);
        a(i,k)=sqrt(g)*cos(theta(i,k)/2);    
        b(i,k)=sqrt(g)*sin(theta(i,k)/2);    
    end
    
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function  [cm, Lm, am, bm]=cmedia(cent,L,a,b,n,N) % Funcion para calcular las medias directas 

%Datos de entrada: - matriz que almacena los centroides de cada curva
%                  - par (a,b)
%                  - numero de curvas
%                  - numero de áangulos que se toman del intervalo [0,2pi]
%Salida: - centroide medio
%        - par (a,b) medio

cm=zeros(2,1);
am=zeros(N,1);
bm=zeros(N,1);

for i=1:N
    
    am(i)=sum(a(i,:))/n;    %Sumamos por filas                 
    bm(i)=sum(b(i,:))/n;
end

cm(1,1)=sum(cent(1,:))/n;
cm(2,1)=sum(cent(2,:))/n;

Lm=sum(log(L))/n;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [a0m, b0m]=nuevosayb(am,bm,N,h) %Funcion para calcular las rotaciones del par (am,bm)

%Datos de entrada: - par (am,bm)
%                  - numero de angulos que se toman del intervalo [0,2pi]
%                  - tamanio de paso 
%Salida: par rotado (a0m,b0m)

%Para calcular este par debemos tener en cuenta varias cosas. Primero, si
%am y bm son ortogonales, entonces a0m y b0m son am y bm normalizados. En
%caso contrario, obtenemos unos nuevos vectores u1 y u2 aplicando
%Gramn-Schimdt a am y bm y los normalizamos, dando como resultado el par (ag,bg).
%Este nuevo par lo rotamos multiplicandolo por e^(ialpha), donde alpha
%varia entre [0,2pi]. Por ultimo, se calcula el alpha que minimiza la
%distancia d0((am,bm),(ag*cos(alpha), bg*sen(alpha)) y se obtienen el a0m y
%b0m.

a0m=zeros(N,1);
b0m=zeros(N,1);

q=integral(h,am.*bm);
%producto interior <am,bm> en V
  
normam=sqrt(integral(h,am.^2));
normbm=sqrt(integral(h,bm.^2));
       
if abs(q)<=1.e-2
    
    %Como estamos aplicando unn metodo numerico, no podemos  exigir que q
    %sea exactamente cero, sino menor que el minimo error del metodo 
    %Trapezoidal.
    
    a0m=am/normam;
    b0m=bm/normbm;
    
    disp('am y bm ya eran ortogonales.');
    
else   %Si am y bm no son ortogonales
    
    u1=am;
    
    u2= bm - (q/normam^2)*am;
    
    normu1=sqrt(integral(h,u1.^2));
    normu2=sqrt(integral(h,u2.^2));

    ag=u1/normu1;
    bg=u2/normu2;
    
    w=integral(h,ag.*bg);   %Comprobacion de que ag y bg son ortonormales
    
    if abs(w)<1.e-2
        
        f1=integral(h,bm.*ag);
        f3=integral(h,am.*ag);
        f4=integral(h,bm.*bg);
        
        if f1 > 0
            
            alpha = atan(f1/(f3+f4));
            
        end
        
        if f1 < 0
            
            alpha =2*pi + atan(f1/(f3+f4));
           
        end
        
        if cos(alpha) > 0
        
            a0m=ag.*cos(alpha) - bg.*sin(alpha);
            b0m=ag.*sin(alpha) + bg.*cos(alpha);
            
            disp('Hemos calculado los nuevos a0m y b0m');
       
       else   %En caso de que ni am y bm, ni ag y bg sean ortogonales
        
            disp('alpha es un máximo');
        
        end
       
    else 
        disp('ag y bg no son ortonormales. FALLO');
    end
end 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function curvaphi2=curvaporphi2(a0m,b0m,N,h)  %Funcion para calcular phi2(a0m,b0m)

%Datos de entrada: - par rotado (a0m, b0m)
%                  - numero de angulos que se toman del intervalo [0,2pi]
%                  - tamanio de paso
%Salida: - matriz que almacena los valores de la curva phi2(a0m,b0m)

%Como tenemos que calcular la integral del numero complejo a0m + ib0m al
%cuadrado entre 0 y theta, tratamos la parte real y la parte imaginaria por serparado,
%juntando los valores en la matriz de salida al final, y truncando las
%evaluaciones. 

curvaphi2=zeros(2,N);
s=zeros(2,N);

for i=1:N
    
    s(1,i)=a0m(i)^2 -b0m(i)^2;
    s(2,i)=2*a0m(i)*b0m(i);
    
    mt=zeros(2,i);
    
    for k=1:i
        
        mt(:,k)=s(:,k);
        
    end
        
    curvaphi2(1,i)= integral(h,mt(1,:))/2; 
    curvaphi2(2,i)= integral(h,mt(2,:))/2;
    
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function d2=der2(curvaphi2,N,h)  % Funcion para calcular la norma de la derivada de la curva phi2(a0m,b0m)

%Datos de entrada: - curva phi2(a0m,b0m)
%                  - numero de angulos que se toman del intervalo [0,2pi]
%                  - tamanio de paso
%Salida: vector con la norma de la derivada de la curva phi2(a0m,b0m)

%Esta derivada se calcula igual que la derivada de las curvas en la
%funcion [d,theta]

der=zeros(2,N);
d2=zeros(N,1);

for i=3:N-2
    der(:,i)=(-curvaphi2(:,i+2)+8*curvaphi2(:,i+1)-8*curvaphi2(:,i-1)+curvaphi2(:,i-2))/(12*h);
end

der(:,1)=(-curvaphi2(:,3)+8*curvaphi2(:,2)-8*curvaphi2(:,N-1)+curvaphi2(:,N-2))/(12*h);
der(:,2)=(-curvaphi2(:,4)+8*curvaphi2(:,3)-8*curvaphi2(:,1)+curvaphi2(:,N-1))/(12*h);
der(:,N-1)=(-curvaphi2(:,2)+8*curvaphi2(:,N)-8*curvaphi2(:,N-2)+curvaphi2(:,N-3))/(12*h);
der(:,N)=(-curvaphi2(:,3)+8*curvaphi2(:,2)-8*curvaphi2(:,N-1)+curvaphi2(:,N-2))/(12*h);

for i=1:N
    d2(i)=sqrt(der(1,i)^2 + der(2,i)^2);
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Longitud=long(d,h) % Funcion para obtener la longitud de una curva 

%Datos de entrada: - derivada de la curva phi2(a0m,b0m)
%                  - tamanio de paso h
%Salida: - longitud de la curva phi2(a0m,b0m)

Longitud= integral(h,d);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function centro=centroi(curva,Long,der,N,h) % Funcion para calcular el centroide de una curva

%Datos de entrada: - curva phi2(a0m,b0m)
%                  - longitud de la curva phi2(a0m,b0m)
%                  - norma de la derivada de la curva phi2(a0m,b0m)
%                  - numero de angulos que se toman del intervalo [0,2pi]
%                  - tamanio de paso
%Salida: - centroide de la curva phi2(a0m,b0m)

%Este centroide se calcula de la misma manera que los centroides de las
%curvas iniciales en la funcion cent.

centro=zeros(2,1);
suma=zeros(2,N);

for i=1:N
    
    suma(:,i)=curva(:,i)*der(i);
  
end

centro(1)=integral(h,suma(1,:))/Long;
centro(2)=integral(h,suma(2,:))/Long;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function curvamedia=curvanueva(cm,Lm,curvaphi2,cent2,N)  %Funcion para calcular la curva media

%Datos de entrada: - centroide medio
%                  - longitud media
%                  - curva phi2(a0m,b0m)
%                  - centroide de la curva phi2(a0m,b0m)
%                  - numero de angulos que se toman del intervalo [0,2pi]
%                  - numero de curvas
%Salida: -curva meddia del conjunto de curvas iniciales.

%La curva media se calcula con la siguiente ecuacion 
% cmedia(t)= cm + (L(1)*...*L(n))^(1/n)*(curvaphi2(a0m,b0m)(t) - cen2)
cm=cm';
cent2=cent2';
curvamedia=zeros(N,2);
curvaphi2=curvaphi2';

for i=1:N
    
    curvamedia(i,:)=cm(1,:) + exp(Lm)*(curvaphi2(i,:) - cent2(1,:));
    
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function I=integral(h,f)  %Funcion Regla de Trapezoidal Compuesta

%Datos de entrada: - tamanio de paso
%                  - funcion de la cual se evalua la integral entre 0 y 2pi
%Salida: - valor de la integral

suma=sum(f);

I=(suma - f(1)/2 - f(end)/2)*h;
 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d3=der3(curvamedia,N,h)  % Funcion para calcular la norma de la derivada de la curva phi2(a0m,b0m)

%Datos de entrada: - curva phi2(a0m,b0m)
%                  - numero de angulos que se toman del intervalo [0,2pi]
%                  - tamanio de paso
%Salida: vector con la norma de la derivada de la curva phi2(a0m,b0m)

%Esta derivada se calcula igual que la derivada de las curvas en la
%funcion [d,theta]

der=zeros(N,2);
d3=zeros(N,1);

for i=3:N-2
    der(i,:)=(-curvamedia(i+2,:)+8*curvamedia(i+1,:)-8*curvamedia(i-1,:)+curvamedia(i-2,:))/(12*h);
end

der(1,:)=(-curvamedia(3,:)+8*curvamedia(2,:)-8*curvamedia(N-1,:)+curvamedia(N-2,:))/(12*h);
der(2,:)=(-curvamedia(4,:)+8*curvamedia(3,:)-8*curvamedia(1,:)+curvamedia(N-1,:))/(12*h);
der(N-1,:)=(-curvamedia(2,:)+8*curvamedia(N,:)-8*curvamedia(N-2,:)+curvamedia(N-3,:))/(12*h);
der(N,:)=(-curvamedia(3,:)+8*curvamedia(2,:)-8*curvamedia(N-1,:)+curvamedia(N-2,:))/(12*h);

for i=1:N
    d3(i)=sqrt(der(i,1)^2 + der(i,2)^2);
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%