#!/usr/bin/octave
# Conjunto de utilidades para explorar distintas redes cristalinas
# Autor: Mariano Marziali Bermúdez
# Versión: 27/3/2019

printf("Conjunto de utilidades para explorar distintas redes cristalinas\n");
printf("Versión: 20/3/2019 (Mariano Marziali Bermúdez)\n");


# Calcula el volumen de la celda primitiva dada por R
function V = volumen3d(R)
    V = abs(det(R));
endfunction

# Calcula una base del espacio recíproco de la red generada por R
function Q = reciproco3d(R)
    Q = (R \ (2 * pi * eye(3)))';
endfunction

# Grafica la red generada por 'R' con base 'S' en la figura 'fig'. 
# 'dim' es un vector que especifica el tamaño de la red [x,y,z]
function h = verRed3d(R, S, dim, fig)
    figure(fig);
    clf; hold on

    L = ceil(2 * dim ./ max(abs(R)));
    n1 = [-L(1):L(1)];
    n2 = [-L(2):L(2)];
    n3 = [-L(3):L(3)];
    [nn1,nn2,nn3] = meshgrid(n1,n2,n3);
    RR = [vec(nn1),vec(nn2),vec(nn3)] * R;
    usar = RR(:,1) <= dim(1) & RR(:,1) >= 0 ...
         & RR(:,2) <= dim(2) & RR(:,2) >= 0 ...
         & RR(:,3) <= dim(3) & RR(:,3) >= 0;
    RR = RR(usar,:);

    nbase = rows(S);
    if (nbase > 1)
        cbase = rainbow(nbase);
        for b = 1:nbase
            h = scatter3(RR(:,1) + S(b,1), RR(:,2) + S(b,2), ...
                         RR(:,3) + S(b,3), 20, cbase(b,:), '.');
        endfor
    else
        h = scatter3(RR(:,1),RR(:,2),RR(:,3), 20, [0,0,0], '.');
    endif
    axis('equal');
endfunction

function tabla = vecinos(R, S, dim)
    L = ceil(2 * dim ./ max(abs(R)));
    n1 = [-L(1):L(1)];
    n2 = [-L(2):L(2)];
    n3 = [-L(3):L(3)];
    [nn1,nn2,nn3] = meshgrid(n1,n2,n3);
    RR = [vec(nn1),vec(nn2),vec(nn3)] * R;
    
    nprim = rows(RR);
    nbase = rows(S);
    if (nbase > 1)
        RR = repmat(RR,nbase,1);
        for b = 1:nbase
            RR((1+(b-1)*nprim):(b*nprim),:) += repmat(S(b,:),nprim,1);
        endfor
    endif
    d = sqrt(RR(:,1).^2 + RR(:,2).^2 + RR(:,3).^2);
    dmin = min(d(d > 0));
    d =  round(1e5 * d/dmin) * dmin/1e5;
    nnd = unique(sort(d(d < max(dim))));
    nn = zeros(size(nnd));
    for k = 1:length(nn)
        nn(k) = sum(d <= nnd(k));
    endfor
    n_vec = diff(nn);
    d_vec = nnd(2:end);
    tabla = [[1:length(n_vec)]', n_vec, d_vec];
endfunction

function [qplot, Iplot] = xrdPolvo(R,S,Z,dim,fig)
    n = 500;
    Q = reciproco3d(R);
    dim = 2*pi*dim;

    L = ceil(2 * dim ./ max(abs(Q)));
    n1 = [-L(1):L(1)];
    n2 = [-L(2):L(2)];
    n3 = [-L(3):L(3)];
    [nn1,nn2,nn3] = meshgrid(n1,n2,n3);
    QQ = [vec(nn1),vec(nn2),vec(nn3)] * Q;
    usar = abs(QQ(:,1)) <= dim(1) ...
         & abs(QQ(:,2)) <= dim(2) ...
         & abs(QQ(:,3)) <= dim(3) ;
    QQ = QQ(usar,:);

    # calculo factores de estructura
    sq = sum(exp(i * QQ * S') * diag(Z), 2);  
    Iq = abs(sq).^2;
    q = sqrt(QQ(:,1).^2 + QQ(:,2).^2 + QQ(:,3).^2);

    qplot = linspace(0, (n+1)/(n-1) * min(dim), n);
    Iplot = zeros(size(qplot));

    for k = 2:length(qplot)
        Iplot(k) = sum(Iq(q <= qplot(k)));
    endfor
    Iplot = [0,diff(Iplot)];
    Iplot = n * Iplot/sum(Iplot);
    

    figure(fig);
    clf; hold on    
    plot(qplot/(2*pi), Iplot, '-k')
    xlabel('|q|/2pi');
    ylabel('intensidad');
    
endfunction