% *************************************************************************
%
% Monitorizacao 2D de estruturas/ajustamento por minimos quadrados
%

% *************************************************************************
%
% rede_2d_v2.m
%
% *************************************************************************
%
% Este programa segue de perto a formulacao descrita em:
% 
% Concepts of network and deformation analysis, W.F.Caspary
% Monograph 11, School of Surveying, The New South Wales 
% Kensington, N.S.W. Australia
% Edited by J.M.Ruegger
% Second (corrected) impression, August 1988
%
% *************************************************************************
%
% Na monitorizacao da deformacao, a deteccao de movimentos de uma
% estrutura pode ser verificada apenas com o formalismo dos minimos
% quadrados e com a analise estatistica associada pois neste caso a
% magnitude dos movimentos potenciais e' da mesma ordem de grandeza da
% precisao (incerteza, erro) das observacoes
%
% *************************************************************************
%
% O ajustamento por minimos quadrados e' uma ferramenta para efectuar de
% forma objectiva o controlo de qualidade de um conjunto de observacoes
% redundantes utilizando regras matematicas bem definidas. 
% O objectivo deste programa esta' a montante do problema do ajustamento por
% minimos quadrados: para este problema ter solucao, e' necessario que 
% esteja definido um referencial ou datum; pretende-se entao estudar a 
% influencia que a escolha do datum tem na solucao do problema do 
% ajustamento (calculo dos parametros) e na respectiva precisao
%
% *************************************************************************
%
% O objecto ou area em estudo e' representado por um conjunto de pontos que
% sao monumentalizados de forma perene, cuja seleccao depende da topografia
% e da estrutura do objecto. O numero de pontos da rede depende fortemente
% do proprio objecto e da deformacao prevista. Normalmente distinguem-se os
% pontos de referencia e os pontos objecto, sendo as observacoes repetidas
% em epocas diferentes, cujo ajustamento e' efectuado de forma independente.
% Apos a primeira repeticao, quando se obtem dois conjuntos de parametros
% (coordenadas) ajustados, e' possivel comparar as coordenadas de cada ponto
% e calcular o vector deslocamento associado. Como o quociente sinal/ruido
% e' pequeno, isto e', os deslocamentos eventuais sao da mesma ordem de
% grandeza dos erros de observacao, e' necessario um grande cuidado na
% avaliacao da qualidade das observacoes e na seleccao do modelo matematico
% mais adequado: quaisquer erros no modelo, em especial a presenca de erros
% grosseiros ou sistematicos nao detectados nas observacoes, assim como a
% utilizacao de variancias a priori pouco verosimeis conduzem a
% deformacoes aparentes que contaminam o resultado.
%
% Uma vez que a analise da deformacao utilizando observacoes geodesicas e'
% baseada na comparacao de coordenadas estimadas dos pontos da rede em
% epocas diferentes, a qualidade dessas estimacoes e' de importancia
% fundamental, qualidade essa aferida pela matriz de
% variancias-covariancias dos parametros (coordenadas e/ou cotas)
%
% Os parametros das elipses de erro xc,yc,a,b,teta sao escritos no ficheiro
% elipses.txt, que sera posteriormente lido no microstation atraves da
% macro elipses1.bas (convem lembrar que o teta calculado e escrito
% no ficheiro indicado refere-se ao angulo entre o eixo Y e o semieixo
% maior da elipse, contado no sentido horario; o angulo teta que o
% microstation utiliza refere-se ao angulo entre o eixo X e o semieixo
% maior da elipse, no sentido anti-horario, sendo a conversao entre os dois
% valores de teta efectuada pela macro.
% 
% *************************************************************************
%
% Quantos parametros devem ser fixados? Em numero igual ao deficit de
% caracteristica da matriz A. No caso de o numero de parâmetros fixados ser
% inferior ao deficit de caracteristica da matriz A (rede livre), a 
% inversao da matriz normal N=At*P*A não é possível, não existindo solução
% para o problema; no caso de o numero de parâmetros fixados ser superior 
% ao deficit de caracteristica da matriz A (rede constrangida), 
% introduzem-se distorções na rede devido a erros nas coordenadas dos 
% pontos considerados fixos; no caso de o número de parâmetros fixados ser
% igual ao deficit de caracteristica da matriz A (rede de constrangimento 
% minimo),a soluçao é possível e não há qualquer distorção da estrutura 
% interna da rede. Tem-se entao que o estabelecimento de um datum nao deve 
% afectar a geometria da rede, isto e, a posicao relativa dos pontos deve 
% ser definida apenas pelas observacoes (o datum nao deve impor um 
% constrangimento maior do que o estritamente necessario para tornar o 
% problema do ajustamento possivel, de forma a tornar a matriz N regular); 
% assim, sendo d o deficit de  caracteristica da matriz N (e da matriz A), 
% um datum que fixe um numero de parametros < d mantem a matriz N singular 
% enquanto que um datum que fixe um numero de parametros > d esta' a 
% distorcer a rede, ao obrigar as observacoes a adaptarem-se ao datum.
%
% Os parâmetros que são fixados para tornar o ajustamento possível
% podem ser coordenadas de pontos da rede e/ou valores de distância e/ou 
% rumos entre pontos da rede; outra hipótese, nas redes de traço mínimo, 
% a soluçao e obtida sem fixar as coordenadas de nenhum ponto da rede (o 
% datum corresponde ao centroide da rede e a um rumo e distancia media. 
% Nalgumas aplicações, por exemplo a análise de deformaçao, é vantajoso 
% introduzir o referencial (local) através da selecção dos pontos estáveis
% (em número maior do que o deficit de caracteristica) dentro do conjunto
% de todos os pontos da rede, em vez de considerar como datum todos os pontos 
% da rede; neste caso, o traço que se minimiza é um subconjunto da matriz QX.
% As solucoes de norma minima (total ou parcial) dao origem a elipses de 
% erro de menores dimensoes mas nao sao necessariamente as melhores
% solucoes; isso depende da aplicacao e das circunstancias do projecto
%
% Pode efectuar-se uma mudanca de datum utilizando uma transformacao que
% mantenha a geometria da rede: para que isto aconteca, o numero de
% parametros desta transformacao nao deve exceder d, o que se designa por
% uma transformacao de similaridade
%
% *************************************************************************
%
% Quao "boa" deve ser a aproximacao inicial das coordenadas (parametros)?
% A resposta depende da dimensao da rede, mais especificamente da distancia
% entre pontos. Assim, as coordenadas aproximadas devem refectir a
% geometria da rede; por exemplo, se a distancia entre os pontos A e B for
% igual a 100.0 m, uma incerteza de 0.05 m nas coordenadas(1/2000) e' 
% inaceitavel para redes de alta precisao; no caso de a distancia entre os
% dois pontos ser 10000.00 m, um erro relativo de 1/200000 (mantendo a
% incerteza de 0.05 m) seria aceitavel
%
% *************************************************************************
%
% No ajustamento por minimos quadrados utilizando o metodo parametrico ou
% das equacoes de observacao, as matrizes A (de configuracao do sistema de
% equacoes ou seja, coeficientes dos parametros) e o vector de fecho vao
% sendo actualizadas durante o processo iterativo. 
%
% *************************************************************************
%
% Qual e o objectivo do teste do qui quadrado para a razao das variancias?
% Pretende-se com o teste verificar se os valores das variancias a priori 
% das observacoes sao realistas
%
% No caso de existirem observacoes angulares, de distancia e de rumo, convem
% comecar o ajustamento verificando se as obs. angulares sao aceitaveis;
% para o efeito, corre-se o programa com a opcao 1 apenas com obs.
% angulares, pelo que e necessario fixar 2 pontos, escolhidos por criterios
% geometricos, de forma a o datum ficar simetricamente bem colocado
% relativamente a disposicao dos pontos da rede. 
% Supondo que o ajustamento converge, para alem da informacao relativa aos 
% valores dos parametros e respectiva precisao e de outra informacao 
% estatistica, e' efectuado o teste da razao de varianacias a priori e a 
% posteriori para verificar se estes valores sao estatisticamente 
% semelhantes, caso em que o ajustamento e aceite.  
% Se o ajustamento for rejeitado (caso mais geral), significa
% que o valor utilizado para a precisao das observacoes angulares e' 
% demasiado optimista (provavelmente o que constava do manual do aparelho),
% sendo necessario repetir o ajustamente com um valor para a precisao das 
% observacoes angulares mais alto (isto e, degrada-se a precisao) ate' o 
% teste ser aceite, correspondendo o ultimo valor utilizado para a precisao
% a precisao efectiva das observacoes angulares, que sera utilizado dai 
% para a frente (caso o valor inicial da precisao das observacoes angulares
% fosse  aceite, este valor ia sendo diminuido (para um valor mais preciso)
% ate' o teste ser rejeitado, adoptando-se como precisao das observacoes 
% angulares o ultimo valor que tinha passado no teste.
%
% Resumindo:
%         sigma02=1 ->
% ----->  sigma_diraz -> 
% |       sigma_diraz**2 ->
% |       q_diraz=sigma_diraz**2/sigma02 ->
% |       p_diraz=1/q_diraz ->
% |       so2 ->
% |<----  teste da razao de variancias
%
% De seguida juntam-se as observacoes de distancia e/ou de rumo, cujo peso
% e calculado a partir de
%
% sigma02=p_direccao*sigma_diraz**2=
%        =p_distancia*sigma_dist**2=
%        =p_rumo*sigma_rumo**2
%        => p_distancia=sigma02/sigma_dist**2
%
% e, ainda na mesma opcao, ajusta-se a rede
%
% *************************************************************************
%
% A teoria classica de erros distingue erros aleatorios, erros
% sistematicos e erros grosseiros. 
% Os erros aleatorios sao inevitaveis, imprevisiveis e seguem determinadas 
% regras estatisticas, correspondendo geralmente a pequenas diferencas 
% entre as observacoes e o respectivo valor medio (na estimacao geodesica e
% nos testes de hipoteses, considera-se que seguem uma distribuicao normal)
% Os erros sistematicos definem-se como a diferenca entre o modelo
% funcional e a realidade, podendo em teoria ser eliminados refinando o
% modelo matematico, seja atraves do modelo funcional (introduzindo
% parametros adicionais ou removendo a respectiva influencia reduzindo
% convenientemente as observacoes), seja atraves do modelo estocastico
% (usando correlacoes a priori entre as observacoes ou parametros
% aleatorios).
% Os erros grosseiros resultam do funcionamento defeituoso do equipamento
% ou de falta de cuidado do operador, supondo-se que podem ser evitados
% Esta classificacao e' importante para a utilizacao dos modelos de
% ajustamento mas nao serve para o desenvolvimento de estrategias de
% deteccao de erros grosseiros pois nao e' possivel separar os erros de
% acordo com essa classificacao em aplicacoes praticas.
% O processo de estimacao fornece residuos, que sao uma mistura de todo o
% tipo de erros: define-se outlier como um residuo que, de acordo com
% alguma regra, e' demasiado grande, podendo admitir-se que os outliers
% correspondem a erros grosseiros.
%
% No caso do teste global (teste da razao de variancias) indicar 
% discrepancias entre o modelo e a realidade (rejeicao) e na hipotese
% destas discrepancias estarem relacionadas com o modelo funcional, estas
% podem ser tomadas em consideracao somando ao vector das observacoes um
% vector D: l_corrigido=l+D, obtendo-se apos o ajustamento:
% 
% x_corrigido_ajust=N*A_t*P*(l+D)=x_ajust+N*A_t*P*D
% v_corrigido=v-Q_v*P*D
% v_corrigido_t*P*v_corrigido=v_t*P*v-D_t*P*Q_v*P*D-2*D_t*P*Q_v*P*v_corrigido
% sigma02=E(sigma_posteriori_2)-E(D_t*P*Q_v*P*D)/f
%
% o que significa que, no caso de um modelo incorrecto,  
% E(sigma_posteriori_2) e sempre maior do que sigma02
%
% Metodo dinamarques
% Admite-se que residuos grandes indicam observacoes menos exactas e
% vice-versa; assim, apos o ajustamento, os pesos a priori sao substituidos
% por outros, funcao dos residuos; efectua-se entao novo ajustamento, que
% conduz a novos residuos, a partir dos quais se obtem novos pesos,
% repetindo-se este procedimento ate ser obtida convergencia, tipicamente
% 5-10 iteracoes; e um metodo puramente heuristico, nao sendo assumidas
% quaisquer hipoteses sobre as propriedades estocasticas das observacoes e
% portanto nao sendo efectiados testes estatisticos.
%
% Formato do ficheiro de input para o programa rede_2d:
%
% 1. cada linha contem 2 ou mais campos com comprimento igual a 10 caracteres,
% justificados a direita, separados por espacos em branco
% 2. o primeiro campo contem um comando que, consoante o caso tem 1 ou mais
% quantificadores
% 3. os comandos sao escritos em letras maiusculas
% 4. nao utilizar TABS para separar campos ou dentro de cada campo
% 5. identificar os pontos sempre da mesma forma (por exemplo P1, p1 e P 1
% sao considerados pontos diferentes
%
% Comandos e respectivos qualificadores:
%
% UANG: unidade angular
%       GMS: graus sexagesimais (grau minuto segundo)
%       GDEC: graus decimais
%       GRD: grados
%
% MAXIT: numero maximo de iteracoes no processo iterativo
%
% CONV: criterio de convergencia (diferenca entre iteracoes sucessivas)
%
% NSIGN: nivel de significancia do teste de hipoteses
% (probabilidade de um erro de primeira especie ou probabilidade de
% rejeitar a hipotese nula mesmo sendo esta correcta; a probabilidade de um
% erro de segunda especie ou probabilidade de aceitar a hipotese nula mesmo
% sendo esta incorrecta nao e geralmente calculado). A rejeicao da hipotese
% nula e admitido como uma indicacao de que alguma coisa esta errada no
% ajustamento(modelo funcional ou modelo estocastico)
%
% PPM: componente dependente da distancia, em partes por milhao, da
% precisao linear de um distanciometro
%
% ESC: modulo da escala para as elipses de erro (ex: 1000 no caso de 1/1000)
%
% CONF: factor de confianca das elipses (para as elipses de erros CONF=1)
% 
% MP: coordenadas (X,Y) de um ponto; sigmaX, sigmaY
%
% ANGAZ: observacao de angulo azimutal entre 2 direccoes: em, de, para,
% ang, sigma, medido no sentido retrogado de para
% no caso de UANG=GMS ou UANG=GDEC, o sigma vem em segundos
% no caso de UANG=GRD, o sigma vem em mgon (1 grado=1 gon, 1000 mgon=1 gon)
%
% DIRAZ: observacao de direccao azimutal entre 2 pontos: em, para dir,
% sigma
% no caso de UANG=GMS ou UANG=GDEC, o sigma vem em segundos
% no caso de UANG=GRD, o sigma vem em mgon (1 grado=1 gon, 1000 mgon=1 gon)
%
% DISTH: observacao de distancia horizontal entre 2 pontos: em, para, dist,
% sigma (DISTH vem em m e sigma em mm)
%
% RUMO: observacao de rumo entre 2 pontos: em, para, rumo, sigma
%
% exemplo de ficheiro de input (retirar %): rede_2d_montsalvens_1.txt
%
%      UANG       GRD
%     MAXIT        10
%     CONV     0.002
%     NSIGN      0.05
%       PPM         0
%       ESC     10000
%      CONF         1
%        MP        P1  100.1030  100.0110
%        MP        P2  111.6010  109.0030
%        MP        P3  122.1810  144.0130
%        MP        P4  116.6920  168.0140
%        MP        P6   87.6610  134.1990
%        MP        P7   88.8550  106.2100
%        MP        P9  129.5510  161.8670
%        MP       P10  102.4480   90.1670
%        MP       P11  126.6760   96.8140
%        MP       P12  143.9770  115.7720
%        MP       P13  145.6870  140.4290
%        MP       P14  133.6100  163.0790
%     DIRAZ        P1        P2   0.00000       0.3
%     DIRAZ        P1       P12  20.29698       0.3
%     DIRAZ        P1       P11  49.87515       0.3
%     DIRAZ        P1       P10 127.36456       0.3
%     DIRAZ        P1        P7 274.31913       0.3
%     DIRAZ        P1        P6 320.03257       0.3
%     DIRAZ        P1        P4 357.48467       0.3
%     DIRAZ        P1        P9 370.53894       0.3
%     DIRAZ        P1        P3 371.85831       0.3
%     DIRAZ        P1       P14 373.34232       0.3
%     DIRAZ        P1       P13 396.07187       0.3
%     DIRAZ        P2        P1   0.00000       0.3
%     DIRAZ        P2        P7  34.47395       0.3
%     DIRAZ        P2        P6  93.87969       0.3
%     DIRAZ        P2        P4 147.73200       0.3
%     DIRAZ        P2        P3 160.93590       0.3
%     DIRAZ        P2        P9 163.09130       0.3
%     DIRAZ        P2       P14 166.85981       0.3
%     DIRAZ        P2       P13 194.83639       0.3
%     DIRAZ        P2       P12 229.13432       0.3
%     DIRAZ        P2       P11 285.53857       0.3
%     DIRAZ        P2       P10 371.04924       0.3
%     DIRAZ        P3        P1   0.00000       0.3
%     DIRAZ        P3        P7  16.39492       0.3
%     DIRAZ        P3        P6  52.76034       0.3
%     DIRAZ        P3        P4 156.08095       0.3
%     DIRAZ        P3        P9 195.31749       0.3
%     DIRAZ        P3       P14 204.77238       0.3
%     DIRAZ        P3       P13 280.02676       0.3
%     DIRAZ        P3       P12 328.54962       0.3
%     DIRAZ        P3       P11 364.34953       0.3
%     DIRAZ        P3        P2 389.07668       0.3
%     DIRAZ        P3       P10 392.75631       0.3
%     DIRAZ        P4        P1   0.00000       0.3
%     DIRAZ        P4        P7  11.70847       0.3
%     DIRAZ        P4        P6  29.93088       0.3
%     DIRAZ        P4       P14 302.83602       0.3
%     DIRAZ        P4        P9 313.15495       0.3
%     DIRAZ        P4       P13 333.18105       0.3
%     DIRAZ        P4       P12 354.12588       0.3
%     DIRAZ        P4        P3 370.45360       0.3
%     DIRAZ        P4       P11 375.89797       0.3
%     DIRAZ        P4        P2 390.24610       0.3
%     DIRAZ        P4       P10 396.28848       0.3
%     DIRAZ        P6        P1   0.00000       0.3
%     DIRAZ        P6        P7  19.50571       0.3
%     DIRAZ        P6        P4 267.38348       0.3
%     DIRAZ        P6        P3 304.58672       0.3
%     DIRAZ        P6        P2 373.84725       0.3
%     DISTH        P1        P2   14.5967       0.3
%     DISTH        P1        P3   49.2302       0.3
%     DISTH        P1        P4   69.9973       0.3
%     DISTH        P2        P3   36.5734       0.3
%     DISTH        P2        P4   59.2302       0.3
%     DISTH        P3        P4   24.6209       0.3
%
% *************************************************************************
%
% Funcoes necessarias para o calculo:
%
% dist_hor_cal(X1,Y1,X2,Y2), que calcula a distancia horizontal entre os 
% pontos (X1,Y1) e (X2,Y2)
%
% dist_hor_lin(X1,Y1,X2,Y2), que calcula os coeficientes da matriz
% de configuracao correspondentes a uma equacao de distancia horizontal
% linearizada entre os pontos (X1,Y1) e (X2,Y2)
%
% ang_az_lin(X1,Y1,X2,Y2,X3,Y3), que calcula os coeficientes da matriz
% de configuracao correspondentes a uma equacao de angulo azimutal 
% linearizada com vertice em (X1,Y1), sendo a primeira pontaria para 
% (X2,Y2) e a segunda pontaria para (X3,Y3), supondo que a luneta roda no
% sentido retrogado
%
% rumo_lin(X1,Y1,X2,Y2), que calcula os coeficientes da matriz
% de configuracao correspondentes a uma equacao de rumo linearizada entre
% os pontos (X1,Y1) e (X2,Y2)
%
% *************************************************************************

clc 

clear all;

%format long;
format short

disp '*********************************************************************'
disp '*                                                                   *'
disp '*                 Ajustamento de uma rede livre 2D                  *'
disp '*                    com constrangimento minimo                     *'
disp '*                                                                   *'
disp '*********************************************************************'

disp '                                                          '

disp '(redes livres ou redes locais, como oposiçao a redes  de adensamento,'
disp 'nas quais o datum já está implicitamente definido através da inclusão'
disp 'de pontos de apoio (fixos) de ordem superior, em numero suficiente;'
disp 'numa rede livre, o datum é definido por um subconjunto dos pontos da'
disp 'da rede)'

disp '    '

% leitura do nome do ficheiro de dados ************************************

ficheiro=input ('Ficheiro de dados: ', 's');

fprintf ('  \n')

% abertura do ficheiro de dados

fid=fopen (ficheiro,'r');

% *************************************************************************
% ciclo de leitura do ficheiro de dados ***********************************
% *************************************************************************

j_mp=1; % indice de pontos MP
j_angaz=1; % indice de observacoes de angulo azimutal
j_disth=1; % indice de observacoes de distancia horizontal
j_rumo=1; % indice de observacoes de rumo
j_diraz=1; % indice de observacoes de direccao azimutal 

% enquanto nao se atingir o fim do ficheiro

while ~feof(fid)

% colocar o conteudo da linha na variavel caracter linha 
% (a variavel linha tem comprimento variavel, conforme a estrutura
% de cada linha, sendo lidos todos os caracteres da linha ate' ao ultimo 
% caracter em branco, nao podendo portanto ser acedidos indices maiores
% do que o comprimento da linha)
    
    linha=fgets(fid);
    
% identificacao do comando no primeiro registo de cada linha e extracao das
% variaveis

% UANG = unidade angular
  
    if (linha(7)==85 & linha(8)==65 & linha(9)==78 & linha(10)==71)
        if (linha(17)==71)
            uang='gdec';
            disp 'Observacoes angulares em graus decimais'
        end
        if (linha(18)==71 & linha(19)==82)
            uang=' grd';
            disp 'Observacoes angulares em grados'
        end
        if (linha(18)==71 & linha(19)==77)
            uang=' gms';
            disp 'Observacoes angulares em graus sexagesimais'
        end
    end
    
% MAXIT = numero maximo de iteracoes
  
    if (linha(6)==77 & linha(7)==65 & linha(8)==88 & linha(9)==73 & linha(10)==84)
        maxit=str2num(linha(1,11:20));
        disp(['Numero maximo de iteracoes no processo iterativo = ',num2str(maxit)]); 
    end
  
% CONV = criterio de convergencia
  
    if (linha(7)==67 & linha(8)==79 & linha(9)==78 & linha(10)==86)
        conv=str2num(linha(1,11:20));  
        disp(['Criterio de convergencia: |diferenca entre iteracoes sucessivas| < ',num2str(conv)]); 
    end
 
% NSIGN = nivel de significancia do teste de razao das variancias
  
    if (linha(6)==78 & linha(7)==83 & linha(8)==73 & linha(9)==71 & linha(10)==78)
        nsign=str2num(linha(1,11:20));
        disp(['Nivel de significancia do teste de hipoteses = ',num2str(nsign)]); 
    end

% PPM = partes por milhao, componente variavel da incerteza na distancia

    if (linha(8)==80 & linha(9)==80 & linha(10)==77)
        ppm=str2num(linha(1,11:20));
        disp(['PPM = ',num2str(ppm)]); 
    end

% ESC = escala para a representacao das elipses de erro

    if (linha(8)==69 & linha(9)==83 & linha(10)==67)
        esc=str2num(linha(1,11:20));
        disp(['ESC = ',num2str(esc)]); 
    end    
    
% CONF = factor de confianca
  
    if (linha(7)==67 & linha(8)==79 & linha(9)==78 & linha(10)==70)
        conf=str2num(linha(1,11:20));  
        disp(['Factor de confianca das elipses: ',num2str(conf)]); 
    end    
    
% MP = coordenadas planimetricas do ponto
  
    if (linha(9)==77 & linha(10)==80)
        nome=linha(1,11:20);
        ponto(j_mp)={nome};
        M(j_mp)=str2num(linha(1,21:30));
        P(j_mp)=str2num(linha(1,31:40));
        if (length(linha)==42) % as coordenadas deste ponto sao parametros
            sigma_M(j_mp)=-1;
            sigma_P(j_mp)=-1;
        else % as coordenadas deste ponto sao fixas (nao sao parametros)
            sigma_M(j_mp)=str2num(linha(1,41:50));
            sigma_P(j_mp)=str2num(linha(1,51:60));
        end
        j_mp=j_mp+1;
    end
  
    n_mp=j_mp-1;  
  
% ANGAZ = angulo azimutal na estacao entre 2 pontos
  
    if (linha(6)==65 & linha(7)==78 & linha(8)==71 & linha(9)==65 & linha(10)==90)
        em=linha(1,11:20);
        de=linha(1,21:30);
        para=linha(1,31:40);
        for j1=1:j_mp-1
            if em==ponto{j1}
                break
            else
                j1=j1+1;
            end
        end
        for j2=1:j_mp-1
            if de==ponto{j2}
                break
            else
                j2=j2+1;
            end
        end
        for j3=1:j_mp-1
            if para==ponto{j3}
                break
            else
                j3=j3+1;
            end
        end
        em_estacao(j_angaz)=j1;
        de_estacao(j_angaz)=j2;
        para_estacao(j_angaz)=j3;
% se uang=gms, converter para grau e decimal de grau        
        if uang==' gms'
            angaz(j_angaz)=str2num(linha(1,41:50))+str2num(linha(1,51:60))/60+str2num(linha(1,61:70))/3600;
            sigma_angaz(j_angaz)=str2num(linha(1,71:80));          
        else
            angaz(j_angaz)=str2num(linha(1,41:50));
            sigma_angaz(j_angaz)=str2num(linha(1,51:60));
        end
        j_angaz=j_angaz+1;
    end
  
    n_angaz=j_angaz-1;
        
% DIRAZ = direccao azimutal entre 2 pontos
  
    if (linha(6)==68 & linha(7)==73 & linha(8)==82 & linha(9)==65 & linha(10)==90)
        em_dir=linha(1,11:20);
        para_dir=linha(1,21:30);
        for j1=1:j_mp-1
            if em_dir==ponto{j1}
                break
            else
                j1=j1+1;
            end
        end
        for j2=1:j_mp-1
            if para_dir==ponto{j2}
                break
            else
                j2=j2+1;
            end
        end
        em_direccao(j_diraz)=j1;
        para_direccao(j_diraz)=j2;
        if uang==' gms'
            diraz(j_diraz)=str2num(linha(1,31:40))+str2num(linha(1,41:50))/60+str2num(linha(1,51:60))/3600;
            sigma_diraz(j_angaz)=str2num(linha(1,61:70));          
        else
            diraz(j_diraz)=str2num(linha(1,31:40));
            sigma_diraz(j_diraz)=str2num(linha(1,41:50));
        end
        j_diraz=j_diraz+1;
    end
  
    n_diraz=j_diraz-1;

% DISTH = distancia horizontal entre 2 pontos
  
    if (linha(6)==68 & linha(7)==73 & linha(8)==83 & linha(9)==84 & linha(10)==72)
        de_dist=linha(1,11:20);
        para_dist=linha(1,21:30);
        for j4=1:j_mp-1
            if de_dist==ponto{j4}
                break
            else
                j4=j4+1;
            end
        end
        for j5=1:j_mp-1
            if para_dist==ponto{j5}
                break
            else
                j5=j5+1;
            end
        end
        de_distancia(j_disth)=j4;
        para_distancia(j_disth)=j5;
        disth(j_disth)=str2num(linha(1,31:40));
        sigma_disth(j_disth)=str2num(linha(1,41:50));
        j_disth=j_disth+1;     
    end  
  
    n_disth=j_disth-1;

% RUMO = rumo entre entre 2 pontos
  
    if (linha(7)==82 & linha(8)==85 & linha(9)==77 & linha(10)==79)
        de_rum=linha(1,11:20);
        para_rum=linha(1,21:30);
        for j4=1:j_mp-1
            if de_rum==ponto{j4}
                break
            else
                j4=j4+1;
            end
        end
        for j5=1:j_mp-1
            if para_rum==ponto{j5}
                break
            else
                j5=j5+1;
            end
        end
        de_rumo(j_rumo)=j4;
        para_rumo(j_rumo)=j5;
        rumo(j_rumo)=str2num(linha(1,31:40));
        sigma_rumo(j_rumo)=str2num(linha(1,41:50));
        j_rumo=j_rumo+1;     
    end  
  
    n_rumo=j_rumo-1;

end

disp '                                                          '
fprintf ('Fim da leitura do ficheiro de dados.\n')
disp '                                                          '

% *************************************************************************
% fim do ciclo de leitura do ficheiro de dados ****************************
% *************************************************************************

% contagem do numero de pontos estacionados no caso de observacoes de
% direccao azimutal pois neste caso ha um parametro adicional (R0) por cada
% estacionamento do aparelho

ne=1;
if (n_diraz>0)
    estacao=em_direccao(1);
    for i=2:n_diraz
        if (estacao~=em_direccao(i))
            ne=ne+1;
            estacao=em_direccao(i);
        end
    end
end

total_obs=n_angaz+n_disth+n_rumo+n_diraz;

fclose (fid);

fprintf ('Numero de pontos da rede: %u \n', n_mp)
fprintf ('Numero de pontos estacionados: %u \n', ne)
fprintf ('Numero de angulos azimutais observados: %u \n', n_angaz)
fprintf ('Numero de distancias horizontais observadas: %u \n', n_disth)
fprintf ('Numero de rumos observados: %u \n', n_rumo)
fprintf ('Numero de direccoes azimutais observadas: %u \n', n_diraz)

disp '                                                          '

% Pode acontecer ter-se em simultaneo n_angaz>0, n_diraz>0, n_rumo>0,
% n_disth>0, embora nao seja usual isto acontecer; 
% Na generalidade dos casos tem-se n_diraz>o e n_disth>0, nao se 
% considerando observacoes de angulos azimutais e de rumo, no primeiro caso
% porque os angulos azimutais sao obtidos por diferenca de observacoes 
% azimutais, fazendo mais sentido utilizar-se directamente as direccoes 
% azimutais (a nao ser que se pretenda utilizar explicitamente angulos 
% azimutais, por exemplo observados segundo o metodo de Scheiber, entrando 
% as varias combinacoes de angulos entre pontos visados numa dada estacao; 
% tem-se assim n_angaz>0 mas n_diraz=0) e no segundo caso
% porque a observacao de rumos implica a utilizacao de um giroscopio, uma
% hipotese pouco comum; no entanto, caso exista a possibilidade de
% utilizacao de um giroscopio (n_rumo>0), a orientacao da rede fica fixada,
% podendo o giroscopio ser utilizado em todos os estacionamentos (neste
% caso n_diraz=0) ou apenas algumas estacoes (neste caso n_diraz>0.

% variancia a priori por unidade de peso (adimensional)
% funciona como um factor de escala arbitrario atraves do qual a inversa da
% matriz de variancias-covariancias das observacoes e' multiplicada para se
% obter a matriz dos pesos das observacoes; uma vantagem fundamental da
% aplicacao do metodo dos minimos quadrados é a estimacao da variancia a
% baseada nas observacoes e na matriz dos pesos das observacoes, que
% se designa por variancia a posteriori; no caso em que as variancia a
% priori e a posteriori sejam estatisticamente iguais, o ajustamento e
% aceite; em geral, o valor da variancia a priori e considerado igual a 1,
% ou seja, considera-se que a matriz peso das observacoes e igual ao 
% inverso da matriz de variancias-covariancias das observacoes; o nome
% completo - variancia a priori por unidade de peso - deriva do facto de
% que se a variancia de uma observacao for sigma_priori_2 entao o peso
% dessa observacao e igual a 1

% transposicao dos vectores de dados de linhas para colunas

ponto=ponto'; % vector que contem o nome dos pontos
M=M'; % vector que contem a coordenada M dos pontos
P=P'; % vector que contem a coordenada P dos pontos
sigma_M=sigma_M';
sigma_P=sigma_P';

nome=num2str(cell2mat(ponto));

nome1 = cellstr(nome);

if (n_angaz>0)
    em_estacao=em_estacao'; % vector que contem o indice do nome do ponto estacao
    de_estacao=de_estacao'; % vector que contem o indice do nome do 1 ponto visado
    para_estacao=para_estacao'; % vector que contem o indice do nome do 2 ponto visado
    angaz=angaz'; % vector que contem o angulo azimutal definido pelas visadas para os pontos 1 e 2
    sigma_angaz=sigma_angaz'; % vector que contem o desvio padrao do angulo azimutal definido pelas visadas para os pontos 1 e 2
end

if (n_diraz>0)
    em_direccao=em_direccao'; % vector que contem o indice do nome do ponto estacao
    para_direccao=para_direccao'; % vector que contem o indice do nome do 2 ponto visado
    diraz=diraz'; % vector que contem o angulo azimutal definido pelas visadas para os pontos 1 e 2
    sigma_diraz=sigma_diraz'; % vector que contem o desvio padrao do angulo azimutal definido pelas visadas para os pontos 1 e 2
end

if (n_disth>0)
    de_distancia=de_distancia';
    para_distancia=para_distancia';
    disth=disth';
    sigma_disth=sigma_disth';
end

if (n_rumo>0)
    de_rumo=de_rumo';
    para_rumo=para_rumo';
    rumo=rumo';
    sigma_rumo=sigma_rumo';    
end

% conversao de graus ou grados para radianos ***************************

if (uang=='gdec' | uang==' gms' & n_angaz>0)
% conversao dos angulos de graus decimais para radianos
    angaz=angaz*pi/180; % rad
% conversao da precisao angular de segundos para radianos
   sigma_angaz=sigma_angaz/3600*pi/180;
end
if (uang==' grd' & n_angaz>0)
% conversao dos angulos de gons para radianos
    angaz=angaz*pi/200; % rad
% conversao da precisao angular de mgons para radianos
    sigma_angaz=sigma_angaz/1000*pi/200;
end

if ((uang=='gdec' | uang==' gms') & n_diraz>0)
% conversao das direccoes de graus decimais para radianos
    diraz=diraz*pi/180; % rad
% conversao da precisao angular de segundos para radianos
    sigma_diraz=sigma_diraz/3600*pi/180;
end
if (uang==' grd' & n_diraz>0)
% conversao das direccoes de gons para radianos
    diraz=diraz*pi/200; % rad
% conversao da precisao angular de mgons para radianos
    sigma_diraz=sigma_diraz/1000*pi/200;
end

if ((uang=='gdec' | uang==' gms') & n_rumo>0)
    rumo=rumo*pi/180; % rad
% conversao da precisao dos rumos de segundos para radianos
    sigma_rumo=sigma_rumo/3600*pi/180;
end
if (uang==' grd' & n_rumo>0)
% conversao dos rumos de gons para radianos
    rumo=rumo*pi/200; % rad
% conversao da precisao dos rumos de mgons para radianos
    sigma_rumo=sigma_rumo/1000*pi/200;
end

% conversao de unidades lineares para metros ******************************
 
% conversao da precisao das distancias de mm para m
if (n_disth>0)
    sigma_disth=sigma_disth/1000;
end

% vector dos parametros (m)

for j=1:n_mp
    X(2*j-1,1)=M(j);
    X(2*j,1)=P(j);
    j=j+1;
end

X;

% *************************************************************************
% construcao da matriz A, que relaciona as observacoes e os parametros, 
% utilizando as coordenadas aproximadas e as equacoes linearizadas
% *************************************************************************

if (n_diraz>0)
    A=zeros(total_obs,2*n_mp+ne);
else
    A=zeros(total_obs,2*n_mp);
end

% angaz
for j=1:n_angaz;
    func=ang_az_lin(X(2*em_estacao(j)-1),X(2*em_estacao(j)),X(2*de_estacao(j)-1),X(2*de_estacao(j)),X(2*para_estacao(j)-1),X(2*para_estacao(j)));
    A(j,2*em_estacao(j)-1)=func(1,1);
    A(j,2*em_estacao(j))=func(2,1);
    A(j,2*de_estacao(j)-1)=func(3,1);
    A(j,2*de_estacao(j))=func(4,1);
    A(j,2*para_estacao(j)-1)=func(5,1);
    A(j,2*para_estacao(j))=func(6,1);
    j=j+1;
end

% disth
for j=n_angaz+1:n_angaz+n_disth
    func=dist_hor_lin(X(2*de_distancia(j-n_angaz)-1),X(2*de_distancia(j-n_angaz)),X(2*para_distancia(j-n_angaz)-1),X(2*para_distancia(j-n_angaz)));
    A(j,2*de_distancia(j-n_angaz)-1)=func(1,1);
    A(j,2*de_distancia(j-n_angaz))=func(2,1);
    A(j,2*para_distancia(j-n_angaz)-1)=func(3,1);
    A(j,2*para_distancia(j-n_angaz))=func(4,1);
    j=j+1; 
end
    
% rumo
for j=n_angaz+n_disth+1:n_angaz+n_disth+n_rumo
    func=rumo_lin(X(2*de_rumo(j-n_angaz-n_disth)-1),X(2*de_rumo(j-n_angaz-n_disth)),X(2*para_rumo(j-n_angaz-n_disth)-1),X(2*para_rumo(j-n_angaz-n_disth)));
    A(j,2*de_rumo(j-n_angaz-n_disth)-1)=func(1,1);
    A(j,2*de_rumo(j-n_angaz-n_disth))=func(2,1);
    A(j,2*para_rumo(j-n_angaz-n_disth)-1)=func(3,1);
    A(j,2*para_rumo(j-n_angaz-n_disth))=func(4,1);
    j=j+1;
end

% diraz   
for j=n_angaz+n_disth+n_rumo+1:total_obs
    func=direccao_lin(X(2*em_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*em_direccao(j-n_angaz-n_disth-n_rumo)),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)));
    A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(1,1);
    A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo))=func(2,1);
    A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(3,1);
    A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo))=func(4,1);
    j=j+1;
end
% colunas da matriz A que contem o coeficiente do R0 em cada estacao
k=1;
j1=1;
for j=2:n_diraz
    if (em_direccao(j) == em_direccao(j-1))
        A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1;
    else
        A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1;
        k=k+1;
    end 
    j=j+1;
end
if (n_diraz>0)
    A(n_angaz+n_disth+n_rumo+n_diraz,2*n_mp+k)=1;
end

A;

% *************************************************************************
% construcao do vector das observacoes
% *************************************************************************

l=zeros(total_obs,1);

for i=1:n_angaz;
    l(i,1)=angaz(i,1);
    i=i+1;
end

for i=n_angaz+1:n_angaz+n_disth
    l(i,1)=disth(i-n_angaz,1);
    i=i+1;
end

for i=n_angaz+n_disth+1:n_angaz+n_disth+n_rumo;
    l(i,1)=rumo(i-n_angaz-n_disth,1);
    i=i+1;
end

for i=n_angaz+n_disth+n_rumo+1:n_angaz+n_disth+n_rumo+n_diraz;
    l(i,1)=diraz(i-n_angaz-n_disth-n_rumo,1);   
    i=i+1;
end

l;

% *************************************************************************
% construcao da matriz de variancias-covariancias das observacoes
% *************************************************************************

sigma_l=zeros(total_obs,total_obs);

%angaz
for j=1:n_angaz;
    sigma_l(j,j)=sigma_angaz(j)^2;
    j=j+1;
end
    
%dist_h
for j=n_angaz+1:n_angaz+n_disth;
    sigma_l(j,j)=sigma_disth(j-n_angaz)^2;
    j=j+1;
end
    
%rumo
for j=n_angaz+n_disth+1:n_angaz+n_disth+n_rumo;
    sigma_l(j,j)=sigma_rumo(j-n_angaz-n_disth)^2;
    j=j+1;
end

%diraz
for j=n_angaz+n_disth+n_rumo+1:n_angaz+n_disth+n_rumo+n_diraz;
    sigma_l(j,j)=sigma_diraz(j-n_angaz-n_disth-n_rumo)^2;
    j=j+1;
end

sigma_l;

% alteracao dos valores das variancias das observacoes de direccao
% em caso de aceitacao do teste da razao de variancias, melhorar os valores
% de sigma_l ate o teste ser rejeitado, sendo o ultimo valor para o qual  o
% teste foi aceite a precisao efectiva das observacoes; em caso de rejeicao
% do teste da razao de variancias, degradar os valores da sigma_l ate o
% teste ser aceite, sendo este valor a precisao efectiva das observacoes.
% caso existam observacoes de direccao e de distancia, apos ter determinado
% o valor de sigma_posteriori_2, calcular os pesos das observacoes de
% distancia utilizando a expressao
% sigma_posteriori_2=
% p_direccao*sigma_2_direccao=
% p_distancia*sigma_2_distancia=
% p_rumo*sigma_2_rumo
%
% variancia de referencia a priori: ADIMENSIONAL
% a seleccao do valor da variancia de referencia a priori nao influencia a
% estimativa do vector dos parametros nem a respectiva matriz de
% variancias-covariancias, podendo estas quantidades ser obtidas assumindo
% qualquer peso relativo entre as observacoes; no entanto, a variancia de
% referencia a posteriori e influenciada pelo valor assumido para a
% variancia a priori: admitindo o valor 1 para a variancia de referencia a
% priori, se a matriz de variancias-covariancias das observacoes for
% adequado e o modelo matematico representar de forma completa o problema,
% a razao entre a variancia a priori e a variancia a posteriori, obtida
% apos o ajustamento, deve ser proxima da unidade (ou seja, as duas
% variancias sao estatisticamente semelhantes); se isto nao acontecer
% (rejeicao do ajustamento), podem ser consideradas 3 hipoteses: 1)erros
% sistematicos nas observacoes 2)inconsistencia entre o modelo matematico
% utilizado e a realidade 3)matriz de variancias-covariancias das
% observacoes inadequada; neste ultimo caso e necessario reformular os
% pesos das observacoes utilizando uma variancia a priori que altere esta
% situacao

% *************************************************************************
sigma_priori_2=1.0;
% *************************************************************************

% lei dos pesos
% sigma_priori_2=p_dir*sigma_dir_2=p_dist*sigma_dist_2=p_ang*sigma_ang_2=p_rumo*sigma_rumo_2

%matriz cofactor das observacoes

Ql=1/sigma_priori_2*sigma_l;

%matriz peso das observacoes

Pesos=inv(Ql);

%Pesos=eye(49)

% *************************************************************************
    
% calculo da matriz normal

N=A'*Pesos*A;
   
[linhas,colunas]=size (N);
   
dimensao=max(linhas,colunas);
    
caracteristica=rank(N);
    
nulidade=dimensao-caracteristica;

str=['Nulidade da matriz A = ',num2str(nulidade)];
disp (str)
disp '                                                          '                          

if (nulidade>0)
  fprintf('A matriz normal N=A´PA e singular (a rede e livre) => e necessario definir o datum.\n\n'); 
end

% calculo do numero de graus de liberdade da rede

if (n_diraz>0)
    gl=total_obs-(2*n_mp+ne-nulidade); % diraz 0> 1 R0 por cada estacao    
else
    gl=total_obs-(2*n_mp-nulidade);
end
        
str=['(nº de graus de liberdade da rede a ajustar com constrangimento minimo = ',num2str(gl),')'];
disp (str)
disp '                                                              '

% *************************************************************************

disp '*************************************************************************************'
disp '*                                                                                   *'
disp '*    1. Datum com variancia nula                                                    *'
disp '*    2. Datum com constrangimentos                                                  *'
disp '*    3. Datum com traco minimo da matriz cofactor dos parametros                    *'
disp '*    4. Datum com minimizacao parcial do traco da matriz cofactor dos parametros    *'
disp '*                                                                                   *'
disp '*************************************************************************************'
disp '                                                                    '

opcao1=input ('           Opcao: ');

disp '                     '

% *************************************************************************
% *                   OPCAO = 1                                           *
% *************************************************************************

if (opcao1 == 1)
    
    clc
    disp '**************************************************************'
    disp '*                  Datum com variancia nula                  *'
    disp '**************************************************************'    
    disp '                                                              '

% *************************************************************************    
%
% angulos azimutais
%
% 1_1**********************************************************************
    if (n_angaz>0 & n_disth==0 & n_rumo==0 & n_diraz==0)
% *************************************************************************

% se nulidade = 4 significa que 4 colunas da matriz A podem ser obtida como
% com binacao linear das restantes colunas

% estabelecimento do datum fixando dois pontos equivale a eliminar as colunas da matriz A correspondente a esses pontos

        for j=1:n_mp;
            if (j<10)
                str=['        ',num2str(j),' - ',char(ponto(j))];
            else
                str=['       ',num2str(j),' - ',char(ponto(j))];        
            end
            disp (str)
            j=j+1;
        end

        disp '                                                          '

        datum1=input ('Indique o primeiro ponto seleccionado para datum: ');

        disp '                                                          '

        for j=1:n_mp;
            if (j ~= datum1)
                if (j<10)
                    str=['        ',num2str(j),' - ',char(ponto(j))];
                else
                    str=['       ',num2str(j),' - ',char(ponto(j))];        
                end
                disp (str)
                j=j+1;
            end
        end
        
        disp '                                                          '
        
        datum2=input ('Indique  o segundo ponto seleccionado para datum: ');

        disp '                                                          '

        c1=2*datum1-1;
        c2=2*datum1;
        c3=2*datum2-1;
        c4=2*datum2;                                                             
        
% inicio do processo iterativo

        for iteracao=1:maxit
    
            %iteracao

% actualizacao da matriz A

            if (iteracao>1)
                % angaz
                for j=1:n_angaz;
                    func=ang_az_lin(X(2*em_estacao(j)-1),X(2*em_estacao(j)),X(2*de_estacao(j)-1),X(2*de_estacao(j)),X(2*para_estacao(j)-1),X(2*para_estacao(j)));
                    A(j,2*em_estacao(j)-1)=func(1,1);
                    A(j,2*em_estacao(j))=func(2,1);
                    A(j,2*de_estacao(j)-1)=func(3,1);
                    A(j,2*de_estacao(j))=func(4,1);
                    A(j,2*para_estacao(j)-1)=func(5,1);
                    A(j,2*para_estacao(j))=func(6,1);
                    j=j+1;
                end
                % disth
                for j=n_angaz+1:n_angaz+n_disth
                    func=dist_hor_lin(X(2*de_distancia(j-n_angaz)-1),X(2*de_distancia(j-n_angaz)),X(2*para_distancia(j-n_angaz)-1),X(2*para_distancia(j-n_angaz)));
                    A(j,2*de_distancia(j-n_angaz)-1)=func(1,1);
                    A(j,2*de_distancia(j-n_angaz))=func(2,1);
                    A(j,2*para_distancia(j-n_angaz)-1)=func(3,1);
                    A(j,2*para_distancia(j-n_angaz))=func(4,1);
                    j=j+1; 
                end
                % rumo
                for j=n_angaz+n_disth+1:n_angaz+n_disth+n_rumo
                    func=rumo_lin(X(2*de_rumo(j-n_angaz-n_disth)-1),X(2*de_rumo(j-n_angaz-n_disth)),X(2*para_rumo(j-n_angaz-n_disth)-1),X(2*para_rumo(j-n_angaz-n_disth)));
                    A(j,2*de_rumo(j-n_angaz-n_disth)-1)=func(1,1);
                    A(j,2*de_rumo(j-n_angaz-n_disth))=func(2,1);
                    A(j,2*para_rumo(j-n_angaz-n_disth)-1)=func(3,1);
                    A(j,2*para_rumo(j-n_angaz-n_disth))=func(4,1);
                    j=j+1;
                end
                % diraz   
                for j=n_angaz+n_disth+n_rumo+1:total_obs
                    func=direccao_lin(X(2*em_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*em_direccao(j-n_angaz-n_disth-n_rumo)),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)));
                    A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(1,1);
                    A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo))=func(2,1);
                    A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(3,1);
                    A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo))=func(4,1);
                    j=j+1;
                end
                % colunas da matriz A que contem o coeficiente do R0 em cada estacao
                k=1;
                j1=1;
                for j=2:n_diraz
                    if (em_direccao(j) == em_direccao(j-1))
                        A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1; 
                    else
                        A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1;
                        k=k+1;
                    end 
                    j=j+1;
                end
                if (n_diraz>0)
                    A(n_angaz+n_disth+n_rumo+n_diraz,2*n_mp+k)=1;
                end
            end
       
% para eliminar uma coluna, utilizar o operador empty []
            
            A1=A;

            corta=[c1 c2 c3 c4];
            
            A1(:,corta)=[];

            A1;

            A
            A1
            pause
    
            A2=A(:,[c1 c2 c3 c4]);
            
            A2;
     
% calculo da matriz L tal que A1*L=A2 (A1_1 e a inversa generalizada de A1)

            A1_1=pinv(A1);

            L=A1_1*A2;
            
% separacao do vector X em X1 e X2, X2 sendo o datum

            X1=X;

            X1T=transpose(X1);

            X1T(:,corta)=[];
            
            X1=transpose(X1T);
            
            XT=transpose(X);
            
            X2T=[XT(:,[c1 c2 c3 c4])];
            
            X2=transpose (X2T);
    
% calculo dos angulos utilizando as coordenadas aproximadas

            for j=1:n_angaz
                angaz_c(j)=ang_az_cal(X(2*em_estacao(j)-1),X(2*em_estacao(j)),X(2*de_estacao(j)-1),X(2*de_estacao(j)),X(2*para_estacao(j)-1),X(2*para_estacao(j)));
                j=j+1;
            end
    
            angaz_c=angaz_c';
    
% calculo do vector fecho W utilizando os angulos observados e os angulos calculados

            for j=1:n_angaz
                W(j,1)=angaz_c(j)-angaz(j);
                if (W(j,1)<-pi);
                    W(j,1)=W(j,1)+2*pi;
                end        
                if (W(j,1)>pi);
                    W(j,1)=W(j,1)-2*pi;
                end
                j=j+1;
            end
    
            W;

% calculo do acrescimo as coordenadas aproximadas

            dX=-inv(A1'*Pesos*A1)*A1'*Pesos*W;
        
            dX;
           
% calculo das novas coordenadas aproximadas

            if (max(abs(dX)) > conv)
                j=1;
                for i=1:2*n_mp; 
                    if (i~=c1 & i~=c2 & i~=c3 & i~=c4)
                        X(i,1)=X1(j,1)+dX(j,1);
                        j=j+1;
                    end
                    i=i+1;
                end
            else
                disp ('O processo iterativo convergiu ...')
                disp '                                                         ')
                break;
            end
            
            iteracao=iteracao+1;
           
            if (iteracao > maxit)
                disp ('O processo iterativo nao convergiu ...')                
                return
            end
            
        end
        
        fprintf('Numero de iterações:  %d\n', iteracao)

        disp ('   ')

        disp 'Vector dos parametros ajustados:'
        for i=1:n_mp; 
            fprintf('  Coordenadas do ponto %s: (%.3f, %.3f)\n', nome1{i}, X(2*i-1), X(2*i));
            i=i+1;
        end

        disp ('   ')
        
        QX1=inv((A1'*Pesos*A1));

        ii=1;
        for i=1:2*n_mp;
            if (i~=c1 & i~=c2 & i~=c3 & i~=c4)
                jj=1;
                for j=1:2*n_mp;
                    if (j~=c1 & j~=c2 & j~=c3 & j~=c4)
                        QX(i,j)=QX1(ii,jj);
                        jj=jj+1;
                    end
                    j=j+1;
                end
                ii=ii+1;
            end
            i=i+1;
        end
        
%        disp 'Matriz cofactor das coordenadas ajustadas:'
        QX;
        
        tracoQX=0;
        for i=1:2*n_mp;
            tracoQX=tracoQX+QX(i,i);
            i=i+1;
        end
        
        disp 'Traco da matriz cofactor das coordenadas ajustadas:';
        tracoQX


        if (gl > 0)

            disp 'Vector dos residuos das observacoes:'
        
            v=A1*dX+W

            v_segundos=v*180/pi*3600

            Qv=Ql-A1*inv(A1'*Pesos*A1)*A1';
        
            l_adj=l+v;
            for i=1:n_angaz
                l_adj=l_adj*180/pi;
            end
        
            l_adj;
        
            Ql_adj=Ql-Qv;
        
            sigma_posteriori_2=v'*Pesos*v/gl
        
            disp 'Teste da razao de variancias:'
            disp '                             '
 
            qui=icdf('chi2',1-nsign,gl);
        
            if (sigma_posteriori_2*gl/sigma_priori_2 < qui)
                disp 'Aceitar o ajustamento'
%           matriz de variancias-covariancias dos parametros9
                sigma_x=sigma_posteriori_2*QX;
%           matriz de variancias-covariancias dos residuos
                sigma_x=sigma_posteriori_2*Qv;
%           matriz de variancias-covariancias das observacoes ajustadas
                sigma_x=sigma_posteriori_2*Ql_adj;                       
            else
                disp 'Rejeitar o ajustamento'
            end

        else

            disp ('    ')
            disp ('A rede não tem redundancia de observações.')

        end

% elipses de erro

        fid=fopen('c:\ajustamento de observações\elipses1.txt','w');

        for i=1:n_mp
            smaior=sqrt((sigma_posteriori_2*QX(2*i-1,2*i-1)+sigma_posteriori_2*QX(2*i,2*i)+sqrt((sigma_posteriori_2*QX(2*i-1,2*i-1)+sigma_posteriori_2*QX(2*i,2*i))^2-4*(sigma_posteriori_2*QX(2*i-1,2*i-1)*sigma_posteriori_2*QX(2*i,2*i)-sigma_posteriori_2*QX(2*i-1,2*i)^2)))/2)*esc;
            smenor=sqrt((sigma_posteriori_2*QX(2*i-1,2*i-1)+sigma_posteriori_2*QX(2*i,2*i)-sqrt((sigma_posteriori_2*QX(2*i-1,2*i-1)+sigma_posteriori_2*QX(2*i,2*i))^2-4*(sigma_posteriori_2*QX(2*i-1,2*i-1)*sigma_posteriori_2*QX(2*i,2*i)-sigma_posteriori_2*QX(2*i-1,2*i)^2)))/2)*esc;
            if (sigma_posteriori_2*QX(2*i-1,2*i-1) == sigma_posteriori_2*QX(2*i,2*i));
                teta=0;
            elseif (sigma_posteriori_2*QX(2*i-1,2*i-1)==0);
                if (sigma_posteriori_2*QX(2*i-1,2*i-1) > sigma_posteriori_2*QX(2*i,2*i));
                    teta=90;
                else
                    teta=0;
                end
            else
                teta=0.5*atan(2*sigma_posteriori_2*QX(2*i-1,2*i)/(sigma_posteriori_2*QX(2*i,2*i)-sigma_posteriori_2*QX(2*i-1,2*i-1)))*180/pi;
            end
            if (teta < 0)
                teta=teta+180;
            end
            fprintf(fid,'%10.4f %10.4f %10.4f %10.4f %10.4f\n',X(2*i-1),X(2*i),smaior,smenor,teta);
            i=i+1;
        end
        fclose (fid);

% o angulo teta refere-se ao angulo entre o eixo Y e a direccao do
% semi-eixo maior da elipse, entre 0 e 180, no sentido horario; o ficheiro
% elipses.txt, cuja estrutura e (xcentro, ycentro,a,b,teta), vai ser lido
% dentro do microstation pela macro elipses1.bas, que constroi as elipses;
% note-se que em microstation os angulos sao contados a partir do eixo X e
% no sentido anti-horario
%
% elipses1.bas
%
%Sub Main
%
%    Dim answer as String
%    Dim M as Double
%    Dim P as Double
%    Dim a as Double
%    Dim b as Double
%    Dim teta as Double
%    Dim point as MbePoint
% 
%    answer = MbeInputBox ("Ficheiro de input:", "c:\topo_aplic\elipses.txt", "Elipses de erro")
%    
%    MbeSendKeyin "co=green"
%    MbeSendKeyin "wt=0"
%    MbeSendKeyin "lv=10"
%    MbeSendKeyin "lc=0"
%    
%    open answer for input as #1
%  
%    
%    do while not eof(1)
%    
%      input #1,M,P,a,b,teta
%      
%      teta=90.0-teta
%      
%      MbeSendCommand "place ellipse center constrained"
%     
%      point.x = M
%      point.y = P
%      
%      MbeSendDataPoint point
%      
%      point.x = M + a*cos(teta/180.0*PI)
%      point.y = P + a*sin(teta/180.0*PI)
%      
%      MbeSendDataPoint point
%      
%      point.x = M + b*cos((teta+90)/180.0*PI)
%      point.y = P + b*sin((teta+90)/180.0*PI)     
%
%      
%      MbeSendDataPoint point
%      
%      MbeSendReset
%      
%    loop
%    
%    MbeSendReset
%  
%    MbeSendKeyin "fit"
%    MbeSendKeyin "selv 1"
%  
%End sub

        %***
        
    end

% *************************************************************************
%
% direccoes azimutais
%
% 1_2**********************************************************************
    if (n_angaz==0 & n_disth==0 & n_rumo==0 & n_diraz>0)
% *************************************************************************        

% se nulidade = 4 significa que 4 colunas da matriz A podem ser obtida como
% combinacao linear das restantes colunas

% estabelecimento do datum fixando dois pontos equivale a eliminar as colunas da matriz A correspondente a esses pontos

        for j=1:n_mp;
            if (j<10)
                str=['        ',num2str(j),' - ',char(ponto(j))];
            else
                str=['       ',num2str(j),' - ',char(ponto(j))];        
            end
            disp (str)
            j=j+1;
        end

        disp '                                                          '

        datum1=input ('Indique o primeiro ponto seleccionado para datum: ');

        disp '                                                          '
 
        for j=1:n_mp;
            if (j ~=datum1)
                if (j<10)
                    str=['        ',num2str(j),' - ',char(ponto(j))];
                else
                    str=['       ',num2str(j),' - ',char(ponto(j))];        
                end
                disp (str)
                j=j+1;
            end
        end

        disp '                                                          '
        
        datum2=input ('Indique  o segundo ponto seleccionado para datum: ');

        disp '                                                          '

        c1=2*datum1-1;
        c2=2*datum1;
        c3=2*datum2-1;
        c4=2*datum2;                                                             
        
% inicio do processo iterativo

        for iteracao=1:maxit
    
            str=['iteracao = ' ,num2str(iteracao)];
            disp (str)

% actualizacao da matriz A

            if (iteracao>1)
                % angaz
                for j=1:n_angaz;
                    func=ang_az_lin(X(2*em_estacao(j)-1),X(2*em_estacao(j)),X(2*de_estacao(j)-1),X(2*de_estacao(j)),X(2*para_estacao(j)-1),X(2*para_estacao(j)));
                    A(j,2*em_estacao(j)-1)=func(1,1);
                    A(j,2*em_estacao(j))=func(2,1);
                    A(j,2*de_estacao(j)-1)=func(3,1);
                    A(j,2*de_estacao(j))=func(4,1);
                    A(j,2*para_estacao(j)-1)=func(5,1);
                    A(j,2*para_estacao(j))=func(6,1);
                    j=j+1;
                end
                % disth
                for j=n_angaz+1:n_angaz+n_disth
                    func=dist_hor_lin(X(2*de_distancia(j-n_angaz)-1),X(2*de_distancia(j-n_angaz)),X(2*para_distancia(j-n_angaz)-1),X(2*para_distancia(j-n_angaz)));
                    A(j,2*de_distancia(j-n_angaz)-1)=func(1,1);
                    A(j,2*de_distancia(j-n_angaz))=func(2,1);
                    A(j,2*para_distancia(j-n_angaz)-1)=func(3,1);
                    A(j,2*para_distancia(j-n_angaz))=func(4,1);
                    j=j+1; 
                end
                % rumo
                for j=n_angaz+n_disth+1:n_angaz+n_disth+n_rumo
                    func=rumo_lin(X(2*de_rumo(j-n_angaz-n_disth)-1),X(2*de_rumo(j-n_angaz-n_disth)),X(2*para_rumo(j-n_angaz-n_disth)-1),X(2*para_rumo(j-n_angaz-n_disth)));
                    A(j,2*de_rumo(j-n_angaz-n_disth)-1)=func(1,1);
                    A(j,2*de_rumo(j-n_angaz-n_disth))=func(2,1);
                    A(j,2*para_rumo(j-n_angaz-n_disth)-1)=func(3,1);
                    A(j,2*para_rumo(j-n_angaz-n_disth))=func(4,1);
                    j=j+1;
                end
                % diraz   
                for j=n_angaz+n_disth+n_rumo+1:total_obs
                    func=direccao_lin(X(2*em_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*em_direccao(j-n_angaz-n_disth-n_rumo)),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)));
                    A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(1,1);
                    A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo))=func(2,1);
                    A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(3,1);
                    A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo))=func(4,1);
                    j=j+1;
                end
                % colunas da matriz A que contem o coeficiente do R0 em cada estacao
                k=1;
                j1=1;
                for j=2:n_diraz
                    if (em_direccao(j) == em_direccao(j-1))
                        A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1; 
                    else
                        A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1;
                        k=k+1;
                    end 
                    j=j+1;
                end
                if (n_diraz>0)
                    A(n_angaz+n_disth+n_rumo+n_diraz,2*n_mp+k)=1;
                end
            end
  
% para eliminar uma coluna, utilizar o operador empty []

            A1=A;
            
            corta=[c1 c2 c3 c4];
            
            A1(:,corta)=[];

            A1;
    
            A2=A(:,[c1 c2 c3 c4]);
            
            A2;
    
% calculo da matriz L tal que A1*L=A2 (A1_1 e a inversa generalizada de A1)

            A1_1=pinv(A1);

            L=A1_1*A2;
            
% separacao do vector X em X1 e X2, X2 sendo o datum

            X1=X;

            X1T=transpose(X1);

            X1T(:,corta)=[];
            
            X1=transpose(X1T);
            
            XT=transpose(X);
            
            X2T=[XT(:,[c1 c2 c3 c4])];
            
            X2=transpose (X2T);

% transformacao de direccoes azimutais em rumos atraves do calculo do R0 em
% cada estacao; para o efeito, e necessario contar o numero de estacoes
% ocupadas e em cada estacao o numero de pontos visados; em  cada estacao
% calcula-se o rumo ponto estacao-ponto visado, tendo-se entao R0=REV-LEV

% calculo do R0 para cada leitura azimutal
    
            for j=1:n_diraz
                R0(j)=rumo(M(em_direccao(j)),P(em_direccao(j)),M(para_direccao(j)),P(para_direccao(j)))-diraz(j);
                if (R0(j)<0)
                    R0(j)=R0(j)+2*pi;
                end
            end

% calculo do R0 medio e respectivo desvio padrao em cada estacao 

            sum_rumo=R0(1);
            sum_rumo2=sum_rumo^2;
            num_pontarias=1;
            j1=1;
            for j=2:n_diraz
                if (em_direccao(j) == em_direccao(j-1))
                    sum_rumo=sum_rumo+R0(j);
                    sum_rumo2=sum_rumo2+R0(j)^2;
                    num_pontarias=num_pontarias+1;
                else
                    j2=j-1;
                    rumo_medio=sum_rumo/num_pontarias;        
                    std_rumo_medio=sqrt((sum_rumo2-sum_rumo^2/(num_pontarias))/(num_pontarias-1));
                    for k=j1:j2 
                        R0_medio(k)=rumo_medio;
                        sigma_R0_medio(k)=std_rumo_medio;
                        k=k+1;
                    end
                    sum_rumo=R0(j);
                    sum_rumo2=sum_rumo^2;
                    num_pontarias=1;
                    j1=j;
                end 
                j=j+1;
            end

            j2=j-1;
            rumo_medio=sum_rumo/num_pontarias;        
            std_rumo_medio=sqrt((sum_rumo2-sum_rumo^2/(num_pontarias))/(num_pontarias-1));  
            for k=j1:j2 
                R0_medio(k)=rumo_medio;
                sigma_R0_medio(k)=std_rumo_medio;
                k=k+1;
            end

            R0=R0';
            R0_medio=R0_medio';
            sigma_R0_medio=sigma_R0_medio';
    
  % calculo das direccoes azimutais "calculadas", dadas como DIRAZ=R-R0

            j1=1;
            for j=2:n_diraz
                if (em_direccao(j) == em_direccao(j-1))
                    diraz_c(j-1)=rumo_cal(X(2*em_direccao(j-1)-1),X(2*em_direccao(j-1)),X(2*para_direccao(j-1)-1),X(2*para_direccao(j-1)))-R0_medio(j1);
                    if (diraz_c(j-1)<0)
                        diraz_c(j-1)=diraz_c(j-1)+2*pi;
                    end
                else
                   diraz_c(j-1)=rumo_cal(X(2*em_direccao(j-1)-1),X(2*em_direccao(j-1)),X(2*para_direccao(j-1)-1),X(2*para_direccao(j-1)))-R0_medio(j1); 
                   if (diraz_c(j-1)<0)
                       diraz_c(j-1)=diraz_c(j-1)+2*pi;
                   end                   
                   j1=j;
                end 
                j=j+1;
            end
            diraz_c(n_diraz)=rumo_cal(X(2*em_direccao(n_diraz)-1),X(2*em_direccao(n_diraz)),X(2*para_direccao(n_diraz)-1),X(2*para_direccao(n_diraz)))-R0_medio(j1); 
            if (diraz_c(n_diraz)<0)
                diraz_c(n_diraz)=diraz_c(n_diraz)+2*pi;
            end     
                    
            diraz_c=diraz_c';
                       
% calculo do vector fecho W utilizando as direccoes observadas e as direccoes calculadas

            for j=1:n_diraz
                W(j,1)=diraz_c(j)-diraz(j);
                if (W(j,1)<-pi);
                    W(j,1)=W(j,1)+2*pi;
                end        
                if (W(j,1)>pi);
                    W(j,1)=W(j,1)-2*pi;
                end
                j=j+1;
            end
    
            W;
             
% calculo do acrescimo as coordenadas aproximadas
            
            dX=inv(A1'*Pesos*A1)*A1'*Pesos*W;
            
% calculo das novas coordenadas aproximadas

            j=1;
            for i=1:2*n_mp; 
                if (i~=c1 & i~=c2 & i~=c3 & i~=c4)
                    X(i,1)=X1(j,1)+dX(j,1);
                    j=j+1;
                end
                i=i+1;
            end
   
            if (max(abs(dX)) < conv)       
                disp ('                  ')
                disp ('O processo iterativo convergiu ...')
                disp ('                                                ')
                break;
            end
            
            iteracao=iteracao+1;
            
            if (iteracao > maxit)
                disp ('                  ')                
                disp ('O processo iterativo nao convergiu ...')                
                stop
            end            
            
        end

% fim do processo iterativo ...
        
% matriz cofactor dos parametros (incluindo o datum, = 0)

        QX1=inv((A1'*Pesos*A1));

        ii=1;
        for i=1:2*n_mp;
            if (i~=c1 & i~=c2 & i~=c3 & i~=c4)
                jj=1;
                for j=1:2*n_mp;
                    if (j~=c1 & j~=c2 & j~=c3 & j~=c4)
                        QX(i,j)=QX1(ii,jj);
                        jj=jj+1;
                    end
                    j=j+1;
                end
                ii=ii+1;
            end
            i=i+1;
        end
 
% traco da matriz cofactor dos parametros

        tracoQX=0;
        for i=1:2*n_mp;
            tracoQX=tracoQX+QX(i,i);
            i=i+1;
        end
        
        str=['Traco da matriz cofactor das coordenadas ajustadas = ' ,num2str(tracoQX)];
        disp (str)
        disp '         '

% vector dos residuos das observacoes

        if (gl > 0)

            v=A1*dX+W;

            disp 'Vector dos residuos das observacoes:';
            disp '    '
            formatSpec20='%10s%10s   %10.4f mgon';
            for j=1:n_diraz
                str=sprintf(formatSpec20,char(ponto(em_direccao(j))),char(ponto(para_direccao(j))),v(j)/pi*200*1000);
                disp (str)
                j=j+1;
            end
        else
            disp ('    ')
            disp ('A rede não tem redundancia de observações.')
        end

        disp '   '
        
        pause
        
% matriz cofactor dos residuos das observacoes

        Qv=Ql-A1*inv(A1'*Pesos*A1)*A1';

% vector das observacoes ajustadas

        l_adj=l+v;
        
% matriz cofactor das observacoes ajustadas

        Ql_adj=Ql-Qv; 
           
        str=['Variancia a priori = ' ,num2str(sigma_priori_2)];
        disp (str)
        disp '         '

% variancia a posteriori (adimensional):

        sigma_posteriori_2=v'*Pesos*v/gl;

        str=['Variancia a posteriori = ' ,num2str(sigma_posteriori_2)];
        disp (str)
        disp '         '
        
        qui=icdf('chi2',1-nsign,gl);
                   
        disp 'Teste da razao de variancias:'
        disp '                             ' 
        
        if (sigma_posteriori_2*gl/sigma_priori_2 < qui)
            disp '    Aceitar o ajustamento'  
        else
            disp '    Rejeitar o ajustamento'
        end
 
        disp '    '
        
        prompt = 'Quer continuar o ajustamento (s/n)? [s]: ';
        str = input(prompt,'s');
        if isempty(str)
            str = 's';
        end
        
        if (str ~= 's')
            %break
        end
        
        disp '    '

% vector dos parametros
        
        formatSpec1='%10s%10.4f  *%10.4f  *';
        formatSpec2='%10s%10.4f   %10.4f';

        disp 'Coordenadas ajustadas dos pontos da rede:'
        disp ' '
        for i=1:n_mp
            if (2*i-1==c1 | 2*i==c2 | 2*i-1==c3 | 2*i==c4)
                str=sprintf(formatSpec1,char(ponto(i)),(X(2*i-1)),(X(2*i)));                
            else
                str=sprintf(formatSpec2,char(ponto(i)),(X(2*i-1)),(X(2*i)));
            end
            disp (str)
            i=i+1;
        end
        disp '       '
        
        pause
        
        formatSpec3='%11.4g  ';
        disp 'Matriz variancias-covariancias das coordenadas ajustadas:'
        disp ' '
        for i=1:2*n_mp
            str=sprintf(formatSpec3,sigma_posteriori_2*QX(i,:));                
            disp (str)
            i=i+1;
        end
        disp '       '            

        traco_sigmaX=0;
        for i=1:2*n_mp;
            traco_sigmaX=traco_sigmaX+sigma_posteriori_2*QX(i,i);
            i=i+1;
        end
        
        str=['Traco da matriz de variancias-covariancias das coordenadas ajustadas = ' ,num2str(traco_sigmaX),' (medida global de exactidao)'];
        disp (str)
        disp '         '
        
        pause
        
% matriz de variancias-covariancias dos residuos
 
        formatSpec3='%11.4g  ';
        disp 'Matriz variancia-covariancia dos residuos:'
        disp ' '
        for i=1:total_obs
            str=sprintf(formatSpec3,sigma_posteriori_2*Qv(i,:));                
            disp (str)
            i=i+1;
        end
        disp '       '        

        pause

% vector das observacoes ajustadas

        disp 'Vector das observacoes ajustadas:';
        disp '    '
        formatSpec20='%10s%10s   %10.4f gon';
        for j=1:n_diraz
            str=sprintf(formatSpec20,char(ponto(em_direccao(j))),char(ponto(para_direccao(j))),l_adj(j)/pi*200);
            disp (str)
            j=j+1;
        end     
        disp '   '
        
        pause
        
% matriz de variancias-covariancias das observacoes ajustadas        
        
        formatSpec3='%11.4g  ';
        disp 'Matriz variancia-covariancia das observacoes ajustadas:'
        disp ' '
        for i=1:total_obs
            str=sprintf(formatSpec3,sigma_posteriori_2*Ql_adj(i,:));                
            disp (str)
            i=i+1;
        end
        disp '       '        
             
        pause
        
% vector dos rumos ajustados

        for j=1:n_diraz
            R_a(j)=rumo(X(2*em_direccao(j)-1),X(2*em_direccao(j)),X(2*para_direccao(j)-1),X(2*para_direccao(j)));
            if (R_a(j)<0)
                R_a(j)=R_a(j)+2*pi;
            end
            if (R_a(j)>2*pi)
                R_a(j)=R_a(j)-2*pi;
            end
            R_a(j)=R_a(j)*200/pi;            
        end

        disp 'Vector dos rumos ajustados:';
        disp '    '
        formatSpec20='%10s%10s   %10.4f gon';
        for j=1:n_diraz
            str=sprintf(formatSpec20,char(ponto(em_direccao(j))),char(ponto(para_direccao(j))),R_a(j));
            disp (str)
            j=j+1;
        end     
        disp '   '
        
% escrever os resultados em ficheiro(s)

% elipses de erro

        fid=fopen('elipses.txt','w');

        for i=1:n_mp
            smaior=sqrt((sigma_posteriori_2*QX(2*i-1,2*i-1)+sigma_posteriori_2*QX(2*i,2*i)+sqrt((sigma_posteriori_2*QX(2*i-1,2*i-1)+sigma_posteriori_2*QX(2*i,2*i))^2-4*(sigma_posteriori_2*QX(2*i-1,2*i-1)*sigma_posteriori_2*QX(2*i,2*i)-sigma_posteriori_2*QX(2*i-1,2*i)^2)))/2)*esc;
            smenor=sqrt((sigma_posteriori_2*QX(2*i-1,2*i-1)+sigma_posteriori_2*QX(2*i,2*i)-sqrt((sigma_posteriori_2*QX(2*i-1,2*i-1)+sigma_posteriori_2*QX(2*i,2*i))^2-4*(sigma_posteriori_2*QX(2*i-1,2*i-1)*sigma_posteriori_2*QX(2*i,2*i)-sigma_posteriori_2*QX(2*i-1,2*i)^2)))/2)*esc;
            if (sigma_posteriori_2*QX(2*i-1,2*i-1) == sigma_posteriori_2*QX(2*i,2*i));
                teta=0;
            elseif (sigma_posteriori_2*QX(2*i-1,2*i-1)==0);
                if (sigma_posteriori_2*QX(2*i-1,2*i-1) > sigma_posteriori_2*QX(2*i,2*i));
                    teta=90;
                else
                    teta=0;
                end
            else
                teta=0.5*atan(2*sigma_posteriori_2*QX(2*i-1,2*i)/(sigma_posteriori_2*QX(2*i,2*i)-sigma_posteriori_2*QX(2*i-1,2*i-1)))*180/pi;
            end
            if (teta < 0)
                teta=teta+180;
            end
            fprintf(fid,'%10.4f%10.4f%10.4f%10.4f%10.4f\n',X(2*i-1),X(2*i),smaior,smenor,teta);
            i=i+1;
        end

% o angulo teta refere-se ao angulo entre o eixo Y e a direccao do
% semi-eixo maior da elipse, entre 0 e 180, no sentido horario; o ficheiro
% elipses.txt, cuja estrutura e (xcentro, ycentro,a,b,teta), vai ser lido
% dentro do microstation pela macro elipses1.bas, que constroi as elipses;
% note-se que em microstation os angulos sao contados a partir do eixo X e
% no sentido anti-horario
%
% elipses1.bas
%
%Sub Main
%
%    Dim answer as String
%    Dim M as Double
%    Dim P as Double
%    Dim a as Double
%    Dim b as Double
%    Dim teta as Double
%    Dim point as MbePoint
% 
%    answer = MbeInputBox ("Ficheiro de input:", "c:\topo_aplic\elipses.txt", "Elipses de erro")
%    
%    MbeSendKeyin "co=green"
%    MbeSendKeyin "wt=0"
%    MbeSendKeyin "lv=10"
%    MbeSendKeyin "lc=0"
%    
%    open answer for input as #1
%  
%    
%    do while not eof(1)
%    
%      input #1,M,P,a,b,teta
%      
%      teta=90.0-teta
%      
%      MbeSendCommand "place ellipse center constrained"
%     
%      point.x = M
%      point.y = P
%      
%      MbeSendDataPoint point
%      
%      point.x = M + a*cos(teta/180.0*PI)
%      point.y = P + a*sin(teta/180.0*PI)
%      
%      MbeSendDataPoint point
%      
%      point.x = M + b*cos((teta+90)/180.0*PI)
%      point.y = P + b*sin((teta+90)/180.0*PI)     
%
%      
%      MbeSendDataPoint point
%      
%      MbeSendReset
%      
%    loop
%    
%    MbeSendReset
%  
%    MbeSendKeyin "fit"
%    MbeSendKeyin "selv 1"
%  
%End sub
       
        fclose (fid);
        
 % deteccao de outliers      
 
        prompt = 'Quer analizar residuos para deteccao de outliers  (s/n)? [s]: ';
        str = input(prompt,'s');
        
        disp '     '
        
        if isempty(str)
            str = 's';
        end
        
        if (str ~= 's')
            %break
        else            
            for iteracao1=1:maxit
                
                str=['iteracao = ' ,num2str(iteracao1)];
                disp (str)     
                
                for i=1:total_obs
                    if (abs(v(i)*sqrt(Pesos(i,i)/sigma_posteriori_2))<3)
                        Pesos(i,i)=Pesos(i,i);
                    else
                        Pesos(i,i)=Pesos(i,i)*exp(-abs(vi)*sqrt(Pesos(i,i))/(3*sigma_posteriori_2));
                    end
                    i=i+1;
                end

% actualizacao da matriz A

                if (iteracao1>1)
                    % angaz
                    for j=1:n_angaz;
                        func=ang_az_lin(X(2*em_estacao(j)-1),X(2*em_estacao(j)),X(2*de_estacao(j)-1),X(2*de_estacao(j)),X(2*para_estacao(j)-1),X(2*para_estacao(j)));
                        A(j,2*em_estacao(j)-1)=func(1,1);
                        A(j,2*em_estacao(j))=func(2,1);
                        A(j,2*de_estacao(j)-1)=func(3,1);
                        A(j,2*de_estacao(j))=func(4,1);
                        A(j,2*para_estacao(j)-1)=func(5,1);
                        A(j,2*para_estacao(j))=func(6,1);
                        j=j+1;
                    end
                    % disth
                    for j=n_angaz+1:n_angaz+n_disth
                        func=dist_hor_lin(X(2*de_distancia(j-n_angaz)-1),X(2*de_distancia(j-n_angaz)),X(2*para_distancia(j-n_angaz)-1),X(2*para_distancia(j-n_angaz)));
                        A(j,2*de_distancia(j-n_angaz)-1)=func(1,1);
                        A(j,2*de_distancia(j-n_angaz))=func(2,1);
                        A(j,2*para_distancia(j-n_angaz)-1)=func(3,1);
                        A(j,2*para_distancia(j-n_angaz))=func(4,1);
                        j=j+1; 
                    end
                    % rumo
                    for j=n_angaz+n_disth+1:n_angaz+n_disth+n_rumo
                        func=rumo_lin(X(2*de_rumo(j-n_angaz-n_disth)-1),X(2*de_rumo(j-n_angaz-n_disth)),X(2*para_rumo(j-n_angaz-n_disth)-1),X(2*para_rumo(j-n_angaz-n_disth)));
                        A(j,2*de_rumo(j-n_angaz-n_disth)-1)=func(1,1);
                        A(j,2*de_rumo(j-n_angaz-n_disth))=func(2,1);
                        A(j,2*para_rumo(j-n_angaz-n_disth)-1)=func(3,1);
                        A(j,2*para_rumo(j-n_angaz-n_disth))=func(4,1);
                        j=j+1;
                    end
                    % diraz   
                    for j=n_angaz+n_disth+n_rumo+1:total_obs
                        func=direccao_lin(X(2*em_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*em_direccao(j-n_angaz-n_disth-n_rumo)),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)));
                        A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(1,1);
                        A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo))=func(2,1);
                        A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(3,1);
                        A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo))=func(4,1);
                        j=j+1;
                    end
                    % colunas da matriz A que contem o coeficiente do R0 em cada estacao
                    k=1;
                    j1=1;
                    for j=2:n_diraz
                        if (em_direccao(j) == em_direccao(j-1))
                            A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1; 
                        else
                            A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1;
                            k=k+1;
                        end 
                        j=j+1;
                    end
                    if (n_diraz>0)
                        A(n_angaz+n_disth+n_rumo+n_diraz,2*n_mp+k)=1;
                    end
                end
  
% para eliminar uma coluna, utilizar o operador empty []

                A1=A;
            
                corta=[c1 c2 c3 c4];
            
                A1(:,corta)=[];

                A1;
    
                A2=A(:,[c1 c2 c3 c4]);
            
                A2;
    
% calculo da matriz L tal que A1*L=A2 (A1_1 e a inversa generalizada de A1)

                A1_1=pinv(A1);

                L=A1_1*A2;
            
% separacao do vector X em X1 e X2, X2 sendo o datum

                X1=X;

                X1T=transpose(X1);

                X1T(:,corta)=[];
            
                X1=transpose(X1T);
            
                XT=transpose(X);
            
                X2T=[XT(:,[c1 c2 c3 c4])];
            
                X2=transpose (X2T);

% transformacao de direccoes azimutais em rumos atraves do calculo do R0 em
% cada estacao; para o efeito, e necessario contar o numero de estacoes
% ocupadas e em cada estacao o numero de pontos visados; em  cada estacao
% calcula-se o rumo ponto estacao-ponto visado, tendo-se entao R0=REV-LEV

% calculo do R0 para cada leitura azimutal
    
                for j=1:n_diraz
                    R0(j)=rumo(M(em_direccao(j)),P(em_direccao(j)),M(para_direccao(j)),P(para_direccao(j)))-diraz(j);
                    if (R0(j)<0)
                        R0(j)=R0(j)+2*pi;
                    end
                end

% calculo do R0 medio e respectivo desvio padrao em cada estacao 

                sum_rumo=R0(1);
                sum_rumo2=sum_rumo^2;
                num_pontarias=1;
                j1=1;
                for j=2:n_diraz
                    if (em_direccao(j) == em_direccao(j-1))
                        sum_rumo=sum_rumo+R0(j);
                        sum_rumo2=sum_rumo2+R0(j)^2;
                        num_pontarias=num_pontarias+1;
                    else
                        j2=j-1;
                        rumo_medio=sum_rumo/num_pontarias;        
                        std_rumo_medio=sqrt((sum_rumo2-sum_rumo^2/(num_pontarias))/(num_pontarias-1));
                        for k=j1:j2 
                            R0_medio(k)=rumo_medio;
                            sigma_R0_medio(k)=std_rumo_medio;
                            k=k+1;
                        end
                        sum_rumo=R0(j);
                        sum_rumo2=sum_rumo^2;
                        num_pontarias=1;
                        j1=j;
                    end 
                    j=j+1;
                end

                j2=j-1;
                rumo_medio=sum_rumo/num_pontarias;        
                std_rumo_medio=sqrt((sum_rumo2-sum_rumo^2/(num_pontarias))/(num_pontarias-1));  
                for k=j1:j2 
                    R0_medio(k)=rumo_medio;
                    sigma_R0_medio(k)=std_rumo_medio;
                    k=k+1;
                end

                R0=R0';
                R0_medio=R0_medio';
                sigma_R0_medio=sigma_R0_medio';
    
  % calculo das direccoes azimutais "calculadas", dadas como DIRAZ=R-R0

                j1=1;
                for j=2:n_diraz
                    if (em_direccao(j) == em_direccao(j-1))
                        diraz_c(j-1)=rumo_cal(X(2*em_direccao(j-1)-1),X(2*em_direccao(j-1)),X(2*para_direccao(j-1)-1),X(2*para_direccao(j-1)))-R0_medio(j1);
                        if (diraz_c(j-1)<0)
                            diraz_c(j-1)=diraz_c(j-1)+2*pi;
                        end
                    else
                        diraz_c(j-1)=rumo_cal(X(2*em_direccao(j-1)-1),X(2*em_direccao(j-1)),X(2*para_direccao(j-1)-1),X(2*para_direccao(j-1)))-R0_medio(j1); 
                        if (diraz_c(j-1)<0)
                            diraz_c(j-1)=diraz_c(j-1)+2*pi;
                        end                   
                        j1=j;
                    end 
                    j=j+1;
                end
                diraz_c(n_diraz)=rumo_cal(X(2*em_direccao(n_diraz)-1),X(2*em_direccao(n_diraz)),X(2*para_direccao(n_diraz)-1),X(2*para_direccao(n_diraz)))-R0_medio(j1); 
                if (diraz_c(n_diraz)<0)
                    diraz_c(n_diraz)=diraz_c(n_diraz)+2*pi;
                end     
                    
                diraz_c=diraz_c';
                       
% calculo do vector fecho W utilizando as direccoes observadas e as direccoes calculadas

                for j=1:n_diraz
                    W(j,1)=diraz_c(j)-diraz(j);
                    if (W(j,1)<-pi);
                        W(j,1)=W(j,1)+2*pi;
                    end        
                    if (W(j,1)>pi);
                        W(j,1)=W(j,1)-2*pi;
                    end
                    j=j+1;
                end
    
                W;
             
% calculo do acrescimo as coordenadas aproximadas
            
                dX=inv(A1'*Pesos*A1)*A1'*Pesos*W;
            
% calculo das novas coordenadas aproximadas

                j=1;
                for i=1:2*n_mp; 
                    if (i~=c1 & i~=c2 & i~=c3 & i~=c4)
                        X(i,1)=X1(j,1)+dX(j,1);
                        j=j+1;
                    end
                    i=i+1;
                end
                   
                v=A1*dX+W;
                
                if (max(abs(v)) < conv)       
                    disp ('                  ')
                    disp ('O processo iterativo convergiu ...')
                    disp ('                                                ')
                    break;
                end
            
                iteracao1=iteracao1+1;
            
                if (iteracao1 > maxit)
                    disp ('                  ')                
                    disp ('O processo iterativo nao convergiu ...')                
                    break;
                end            
            
            end
            
            disp 'Vector dos residuos das observacoes:';
            disp '    '
            formatSpec20='%10s%10s   %10.4f mgon';
            for j=1:n_diraz
                str=sprintf(formatSpec20,char(ponto(em_direccao(j))),char(ponto(para_direccao(j))),v(j)/pi*200*1000);
                disp (str)
                j=j+1;
            end     
            disp '   '
        
            pause

        end            
        
 % escrita dos resultados num ficheiro
 
        prompt = 'Quer escrever os resultados num ficheiro  (s/n)? [s]: ';
        str = input(prompt,'s');
        if isempty(str)
            str = 's';
        end
        
        if (str ~= 's')
            %break
        end
        
        disp '    '

 
        
  
  
        disp 'Fim do ajustamento!'
    end

% ************************************************************************* 
% distancias horizontais
% 1_3**********************************************************************
    if (n_angaz==0 & n_disth>0 & n_rumo==0 & n_diraz==0)
% *************************************************************************

% se nulidade = 3 significa que 3 colunas da matriz A podem ser obtida como
% com binacao linear das restantes colunas
       
% estabelecimento do datum fixando um ponto mais uma 2coordenadas de um segundo ponto equivale a eliminar as colunas da matriz A correspondente a esses pontos

        for j=1:n_mp;
            if (j<10)
                str=['        ',num2str(j),' - ',char(ponto(j))];
            else
                str=['       ',num2str(j),' - ',char(ponto(j))];        
            end
            disp (str)
            j=j+1;
        end

        disp '                                                          '

        datum1=input ('Indique o primeiro ponto seleccionado para datum: ');
 
        clc
        
        disp '                                                          '
         
        jj=1;
         for j=1:2*n_mp;
              if (j<10 & mod (j,2)~=0)
                 str=['          ',num2str(j),' - X',char(ponto(jj))];
                 disp (str)
              end
              if (j<10 & mod (j,2)==0)
                 str=['          ',num2str(j),' - Y',char(ponto(jj))];
                 disp (str)
              end
              if (j>9 & mod (j,2)~=0)
                 str=['         ',num2str(j),' - X',char(ponto(jj))];
                 disp (str)
              end
              if (j>9 & mod (j,2)==0)
                 str=['         ',num2str(j),' - Y',char(ponto(jj))];
                 disp (str)
              end
              if (mod(j,2)==0)
                  jj=jj+1;
              end
              j=j+1;
         end
                    
        disp '                                                          '
        
        datum2=input ('Indique uma coordenada do segundo ponto seleccionado para datum: ');

        c1=2*datum1-1;
        c2=2*datum1;
        c3=datum2;

% inicio do processo iterativo

        for iteracao=1:maxit
    
            iteracao

% actualizacao da matriz A

            if (iteracao>1)
                % angaz
                for j=1:n_angaz;
                    func=ang_az_lin(X(2*em_estacao(j)-1),X(2*em_estacao(j)),X(2*de_estacao(j)-1),X(2*de_estacao(j)),X(2*para_estacao(j)-1),X(2*para_estacao(j)));
                    A(j,2*em_estacao(j)-1)=func(1,1);
                    A(j,2*em_estacao(j))=func(2,1);
                    A(j,2*de_estacao(j)-1)=func(3,1);
                    A(j,2*de_estacao(j))=func(4,1);
                    A(j,2*para_estacao(j)-1)=func(5,1);
                    A(j,2*para_estacao(j))=func(6,1);
                    j=j+1;
                end
                % disth
                for j=n_angaz+1:n_angaz+n_disth
                    func=dist_hor_lin(X(2*de_distancia(j-n_angaz)-1),X(2*de_distancia(j-n_angaz)),X(2*para_distancia(j-n_angaz)-1),X(2*para_distancia(j-n_angaz)));
                    A(j,2*de_distancia(j-n_angaz)-1)=func(1,1);
                    A(j,2*de_distancia(j-n_angaz))=func(2,1);
                    A(j,2*para_distancia(j-n_angaz)-1)=func(3,1);
                    A(j,2*para_distancia(j-n_angaz))=func(4,1);
                    j=j+1; 
                end
                % rumo
                for j=n_angaz+n_disth+1:n_angaz+n_disth+n_rumo
                    func=rumo_lin(X(2*de_rumo(j-n_angaz-n_disth)-1),X(2*de_rumo(j-n_angaz-n_disth)),X(2*para_rumo(j-n_angaz-n_disth)-1),X(2*para_rumo(j-n_angaz-n_disth)));
                    A(j,2*de_rumo(j-n_angaz-n_disth)-1)=func(1,1);
                    A(j,2*de_rumo(j-n_angaz-n_disth))=func(2,1);
                    A(j,2*para_rumo(j-n_angaz-n_disth)-1)=func(3,1);
                    A(j,2*para_rumo(j-n_angaz-n_disth))=func(4,1);
                    j=j+1;
                end
                % diraz   
                for j=n_angaz+n_disth+n_rumo+1:total_obs
                    func=direccao_lin(X(2*em_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*em_direccao(j-n_angaz-n_disth-n_rumo)),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)));
                    A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(1,1);
                    A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo))=func(2,1);
                    A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(3,1);
                    A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo))=func(4,1);
                    j=j+1;
                end
                % colunas da matriz A que contem o coeficiente do R0 em cada estacao
                k=1;
                j1=1;
                for j=2:n_diraz
                    if (em_direccao(j) == em_direccao(j-1))
                        A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1; 
                    else
                        A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1;
                        k=k+1;
                    end 
                    j=j+1;
                end
                if (n_diraz>0)
                    A(n_angaz+n_disth+n_rumo+n_diraz,2*n_mp+k)=1;
                end
            end
       
% para eliminar uma coluna, utilizar o operador empty []

            A1=A;

            corta=[c1 c2 c3];
            
            A1(:,corta)=[];

            A1
    
            A2=A(:,[c1 c2 c3]);
            
            A2
     
% calculo da matriz L tal que A1*L=A2 (A1_1 e a inversa generalizada de A1)

            A1_1=pinv(A1);

            L=A1_1*A2;

% separacao do vector X em X1 e X2, X2 sendo o datum

            X1=X;

            X1T=transpose(X1);

            X1T(:,corta)=[];
            
            X1=transpose(X1T);
            
            XT=transpose(X);
            
            X2T=[XT(:,[c1 c2 c3])];
            
            X2=transpose (X2T);
    
% calculo das distancias utilizando as coordenadas aproximadas

            for j=1:n_disth
                dist_c(j)=dist_hor_cal(X(2*de_distancia(j)-1),X(2*de_distancia(j)),X(2*para_distancia(j)-1),X(2*para_distancia(j)));
                j=j+1;
            end
                
% calculo do vector W utilizando as distancias calculadas e as distancias observadas

            for j=n_angaz+1:n_angaz+n_disth
                W(j,1)=dist_c(j-n_angaz)-dist(j-n_angaz);
                j=j+1;
            end
    
            W;
            
% calculo do acrescimo as coordenadas aproximadas

            dX=-inv(A1'*Pesos*A1)*A1'*Pesos*W;
        
            dX;
           
% calculo das novas coordenadas aproximadas

            if (max(abs(dX)) > conv)
                j=1;
                for i=1:2*n_mp; 
                    if (i~=c1 & i~=c2 & i~=c3)
                        X(i,1)=X1(j,1)+dX(j,1);
                        j=j+1;
                    end
                    i=i+1;
                end
            else
                disp ('O processo iterativo convergiu ...')
                break;
            end
            
            iteracao=iteracao+1;
            
            if (iteracao>maxit)
                disp ('O processo iterativo nao convergiu ...')                
                stop
            end              
        end

        disp 'Vector dos parametros ajustados:'
        X
        
        QX1=inv((A1'*Pesos*A1));

        ii=1;
        for i=1:2*n_mp;
            if (i~=c1 & i~=c2 & i~=c3)
                jj=1;
                for j=1:2*n_mp;
                    if (j~=c1 & j~=c2 & j~=c3)
                        QX(i,j)=QX1(ii,jj);
                        jj=jj+1;
                    end
                    j=j+1;
                end
                ii=ii+1;
            end
            i=i+1;
        end
        
%        disp 'Matriz cofactor das coordenadas ajustadas:'
        QX;
        
        tracoQX=0;
        for i=1:2*n_mp;
            tracoQX=tracoQX+QX(i,i);
            i=i+1;
        end
        
        disp 'Traco da matriz cofactor das coordenadas ajustadas:'
        tracoQX
                
        disp 'Vector dos residuos das observacoes:'
        v=A1*dX+W
        
        Qv=Ql-A1*inv(A1'*Pesos*A1)*A1';
        
        l_adj=l+v;
        for i=1:n_angaz
            l_adj=l_adj*180/pi;
        end
        
        l_adj;
        
        Ql_adj=Ql-Qv;
        
        sigma_posteriori_2=v'*Pesos*v/gl;
        
        disp 'Teste da razao de variancias:'
        disp '                             '
 
        qui=icdf('chi2',1-nsign,gl);
        
        if (sigma_posteriori_2*gl/sigma_priori_2 < qui)
            disp 'Aceitar o ajustamento'
        else
            disp 'Rejeitar o ajustamento'
        end
    
    end

% *************************************************************************
% distancias horizontais+direccoes azimutais
% 1_4**********************************************************************
    if (n_angaz==0 & n_disth>0 & n_rumo==0 & n_diraz>0)
% ************************************************************************* 

% se nulidade = 3 significa que 3 colunas da matriz A podem ser obtida como
% com binacao linear das restantes colunas
      
% estabelecimento do datum fixando um ponto mais uma coordenada de um segundo ponto equivale a eliminar as colunas da matriz A correspondente a esses pontos

        for j=1:n_mp;
            if (j<10)
                str=['        ',num2str(j),' - ',char(ponto(j))];
            else
                str=['       ',num2str(j),' - ',char(ponto(j))];        
            end
            disp (str)
            j=j+1;
        end

        disp '                                                          '

        datum1=input ('Indique o primeiro ponto seleccionado para datum: ');
 
        clc
        
        disp '                                                          '
         
        jj=1;
         for j=1:2*n_mp;
              if (j<10 & mod (j,2)~=0 & j~=2*datum1 & j~=(2*datum1-1))
                str=['          ',num2str(j),' - X',char(ponto(jj))];
                 disp (str)
              end
              if (j<10 & mod (j,2)==0 & j~=2*datum1 & j~=(2*datum1-1))
                 str=['          ',num2str(j),' - Y',char(ponto(jj))];
                 disp (str)
              end
              if (j>9 & mod (j,2)~=0 & j~=2*datum1 & j~=(2*datum1-1))
                 str=['         ',num2str(j),' - X',char(ponto(jj))];
                 disp (str)
              end
              if (j>9 & mod (j,2)==0 & j~=2*datum1 & j~=(2*datum1-1))
                 str=['         ',num2str(j),' - Y',char(ponto(jj))];
                 disp (str)
              end
              if (mod(j,2)==0)
                  jj=jj+1;
              end
              j=j+1;
         end
         
        disp '                                                          '
        
        datum2=input ('Indique uma coordenada do segundo ponto seleccionado para datum: ');

        c1=2*datum1-1;
        c2=2*datum1;
        c3=datum2;
               
% inicio do processo iterativo

        for iteracao=1:maxit
    
            iteracao

% actualizacao da matriz A

            if (iteracao>1)
                % angaz
                for j=1:n_angaz;
                    func=ang_az_lin(X(2*em_estacao(j)-1),X(2*em_estacao(j)),X(2*de_estacao(j)-1),X(2*de_estacao(j)),X(2*para_estacao(j)-1),X(2*para_estacao(j)));
                    A(j,2*em_estacao(j)-1)=func(1,1);
                    A(j,2*em_estacao(j))=func(2,1);
                    A(j,2*de_estacao(j)-1)=func(3,1);
                    A(j,2*de_estacao(j))=func(4,1);
                    A(j,2*para_estacao(j)-1)=func(5,1);
                    A(j,2*para_estacao(j))=func(6,1);
                    j=j+1;
                end
                % disth
                for j=n_angaz+1:n_angaz+n_disth
                    func=dist_hor_lin(X(2*de_distancia(j-n_angaz)-1),X(2*de_distancia(j-n_angaz)),X(2*para_distancia(j-n_angaz)-1),X(2*para_distancia(j-n_angaz)));
                    A(j,2*de_distancia(j-n_angaz)-1)=func(1,1);
                    A(j,2*de_distancia(j-n_angaz))=func(2,1);
                    A(j,2*para_distancia(j-n_angaz)-1)=func(3,1);
                    A(j,2*para_distancia(j-n_angaz))=func(4,1);
                    j=j+1; 
                end
                % rumo
                for j=n_angaz+n_disth+1:n_angaz+n_disth+n_rumo
                    func=rumo_lin(X(2*de_rumo(j-n_angaz-n_disth)-1),X(2*de_rumo(j-n_angaz-n_disth)),X(2*para_rumo(j-n_angaz-n_disth)-1),X(2*para_rumo(j-n_angaz-n_disth)));
                    A(j,2*de_rumo(j-n_angaz-n_disth)-1)=func(1,1);
                    A(j,2*de_rumo(j-n_angaz-n_disth))=func(2,1);
                    A(j,2*para_rumo(j-n_angaz-n_disth)-1)=func(3,1);
                    A(j,2*para_rumo(j-n_angaz-n_disth))=func(4,1);
                    j=j+1;
                end
                % diraz   
                for j=n_angaz+n_disth+n_rumo+1:total_obs
                    func=direccao_lin(X(2*em_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*em_direccao(j-n_angaz-n_disth-n_rumo)),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)));
                    A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(1,1);
                    A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo))=func(2,1);
                    A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(3,1);
                    A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo))=func(4,1);
                    j=j+1;
                end
                % colunas da matriz A que contem o coeficiente do R0 em cada estacao
                k=1;
                j1=1;
                for j=2:n_diraz
                    if (em_direccao(j) == em_direccao(j-1))
                        A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1; 
                    else
                        A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1;
                        k=k+1;
                    end 
                    j=j+1;
                end
                if (n_diraz>0)
                    A(n_angaz+n_disth+n_rumo+n_diraz,2*n_mp+k)=1;
                end
            end
       
% para eliminar uma coluna, utilizar o operador empty []

            A1=A;

            corta=[c1 c2 c3];
            
            A1(:,corta)=[];

            A1;
    
            A2=A(:,[c1 c2 c3]);
            
            A2;
     
% calculo da matriz L tal que A1*L=A2 (A1_1 e a inversa generalizada de A1)

            A1_1=pinv(A1);

            L=A1_1*A2;

% separacao do vector X em X1 e X2, X2 sendo o datum

            X1=X;

            X1T=transpose(X1);

            X1T(:,corta)=[];
            
            X1=transpose(X1T)
            
            XT=transpose(X);
            
            X2T=[XT(:,[c1 c2 c3])];
            
            X2=transpose (X2T)
    
% calculo das distancias utilizando as coordenadas aproximadas

            for j=1:n_disth
                dist_c(j)=dist_hor_cal(X(2*de_distancia(j)-1),X(2*de_distancia(j)),X(2*para_distancia(j)-1),X(2*para_distancia(j)));
                j=j+1;
            end
            
            dist_c=dist_c';
              
% transformacao de direccoes azimutais em rumos atraves do calculo do R0 em
% cada estacao; para o efeito, e necessario contar o numero de estacoes
% ocupadas e em cada estacao o numero de pontos visados; em  cada estacao
% calcula-se o rumo ponto estacao-ponto visado, tendo-se entao R0=REV-LEV

% calculo do R0 para cada leitura azimutal
    
            for j=1:n_diraz
                R0(j)=rumo(M(em_direccao(j)),P(em_direccao(j)),M(para_direccao(j)),P(para_direccao(j)))-diraz(j);
                if (R0(j)<0)
                    R0(j)=R0(j)+2*pi;
                end
            end

% calculo do R0 medio e respectivo desvio padrao em cada estacao 

            sum_rumo=R0(1);
            sum_rumo2=sum_rumo^2;
            num_pontarias=1;
            j1=1;
            for j=2:n_diraz
                if (em_direccao(j) == em_direccao(j-1))
                    sum_rumo=sum_rumo+R0(j);
                    sum_rumo2=sum_rumo2+R0(j)^2;
                    num_pontarias=num_pontarias+1;
                else
                    j2=j-1;
                    rumo_medio=sum_rumo/num_pontarias;        
                    std_rumo_medio=sqrt((sum_rumo2-sum_rumo^2/(num_pontarias))/(num_pontarias-1));
                    for k=j1:j2 
                        R0_medio(k)=rumo_medio;
                        sigma_R0_medio(k)=std_rumo_medio;
                        k=k+1;
                    end
                    sum_rumo=R0(j);
                    sum_rumo2=sum_rumo^2;
                    num_pontarias=1;
                    j1=j;
                end 
                j=j+1;
            end

            j2=j-1;
            rumo_medio=sum_rumo/num_pontarias;        
            std_rumo_medio=sqrt((sum_rumo2-sum_rumo^2/(num_pontarias))/(num_pontarias-1));  
            for k=j1:j2 
                R0_medio(k)=rumo_medio;
                sigma_R0_medio(k)=std_rumo_medio;
                k=k+1;
            end

            R0=R0';
            R0_medio=R0_medio';
            sigma_R0_medio=sigma_R0_medio';
    
  % calculo das direccoes azimutais "calculadas", dadas como DIRAZ=R-R0

            j1=1;
            for j=2:n_diraz
                if (em_direccao(j) == em_direccao(j-1))
                    diraz_c(j-1)=rumo_cal(X(2*em_direccao(j-1)-1),X(2*em_direccao(j-1)),X(2*para_direccao(j-1)-1),X(2*para_direccao(j-1)))-R0_medio(j1);
                    if (diraz_c(j-1)<0)
                        diraz_c(j-1)=diraz_c(j-1)+2*pi;
                    end
                else
                   diraz_c(j-1)=rumo_cal(X(2*em_direccao(j-1)-1),X(2*em_direccao(j-1)),X(2*para_direccao(j-1)-1),X(2*para_direccao(j-1)))-R0_medio(j1); 
                   if (diraz_c(j-1)<0)
                       diraz_c(j-1)=diraz_c(j-1)+2*pi;
                   end                   
                   j1=j;
                end 
                j=j+1;
            end
            diraz_c(n_diraz)=rumo_cal(X(2*em_direccao(n_diraz)-1),X(2*em_direccao(n_diraz)),X(2*para_direccao(n_diraz)-1),X(2*para_direccao(n_diraz)))-R0_medio(j1); 
            if (diraz_c(n_diraz)<0)
                diraz_c(n_diraz)=diraz_c(n_diraz)+2*pi;
            end     
                    
            diraz_c=diraz_c';
                       
% calculo do vector W (vector fecho)
           
            for j=n_angaz+1:n_angaz+n_disth
                W(j,1)=dist_c(j-n_angaz)-disth(j-n_angaz);
                j=j+1;
            end
            for j=1:n_diraz
                W(n_disth+j,1)=diraz_c(j)-diraz(j);
                if (W(n_disth+j,1)<-pi);
                    W(n_disth+j,1)=W(n_disth+j,1)+2*pi;
                end        
                if (W(n_disth+j,1)>pi);
                    W(n_disth+j,1)=W(n_disth+j,1)-2*pi;
                end
                j=j+1;
            end
    
            W;
                       
% calculo do acrescimo as coordenadas aproximadas

            dX=inv(A1'*Pesos*A1)*A1'*Pesos*W;
        
            dX
           
% calculo das novas coordenadas aproximadas

            if (max(abs(dX)) > conv)
                j=1;
                for i=1:2*n_mp; 
                    if (i~=c1 & i~=c2 & i~=c3)
                        X(i,1)=X1(j,1)+dX(j,1);
                        j=j+1;
                    end
                    i=i+1;
                end
            else
                disp ('O processo iterativo convergiu ...')
               break;
            end
            
            iteracao=iteracao+1;
            
            if (iteracao>maxit)
                disp ('O processo iterativo nao convergiu ...')                
                stop
            end  
            
        end

        disp 'Vector dos parametros ajustados:'
        X
        
        QX1=inv((A1'*Pesos*A1));

        ii=1;
        for i=1:2*n_mp;
            if (i~=c1 & i~=c2 & i~=c3)
                jj=1;
                for j=1:2*n_mp;
                    if (j~=c1 & j~=c2 & j~=c3)
                        QX(i,j)=QX1(ii,jj);
                        jj=jj+1;
                    end
                    j=j+1;
                end
                ii=ii+1;
            end
            i=i+1;
        end
        
%        disp 'Matriz cofactor das coordenadas ajustadas:'
        QX;
        
        tracoQX=0;
        for i=1:2*n_mp;
            tracoQX=tracoQX+QX(i,i);
            i=i+1;
        end
        
        disp 'Traco da matriz cofactor das coordenadas ajustadas:'
        tracoQX
 
        disp 'Vector dos residuos das observacoes:'        
        v=A1*dX+W
        
        Qv=Ql-A1*inv(A1'*Pesos*A1)*A1';
        
        l_adj=l+v;
        for i=1:n_angaz
            l_adj=l_adj*180/pi;
        end
        
        l_adj;
        
        Ql_adj=Ql-Qv;
        
        sigma_posteriori_2=v'*Pesos*v/gl;
        
        disp 'Teste da razao de variancias:'
        disp '                             '
 
        qui=icdf('chi2',1-nsign,gl);
        
        if (sigma_posteriori_2*gl/sigma_priori_2 < qui)
            disp 'Aceitar o ajustamento'
        else
            disp 'Rejeitar o ajustamento'
        end
    
    end
    
% 1_5*********************************************************************
    if (n_disth > 0 & n_rumo>0)
% ************************************************************************ 
        
% se nulidade = 2 significa que 2 colunas da matriz A podem ser obtida como com binacao linear das restantes colunas

% estabelecimento do datum fixando um ponto equivale a eliminar as colunas da matriz A correspondente a esse ponto

        for j=1:n_mp;
            if (j<10)
                str=['        ',num2str(j),' - ',char(ponto(j))];
            else
                str=['       ',num2str(j),' - ',char(ponto(j))];        
            end
            disp (str)
            j=j+1;
        end

        disp '                                                          '

        datum1=input ('Indique o ponto seleccionado para datum: ');
 
        clc
        
        c1=2*datum1-1;
        c2=2*datum1;

% calculo do numero de graus de liberdade da rede

        if (n_diraz>0)
            gl=total_obs-(2*n_mp-nulidade+ne); % ne=1 R0 por cada estacao
        else
            gl=total_obs-(2*n_mp-nulidade);
        end
        str=['nº de graus de liberdade da rede a ajustar com constrangimento minimo = ',num2str(gl)];
        disp (str)
        disp '                                                          '  
        
% inicio do processo iterativo

        for iteracao=1:maxit
    
            iteracao

% actualizacao da matriz A

            if (iteracao>1)
                % angaz
                for j=1:n_angaz;
                    func=ang_az_lin(X(2*em_estacao(j)-1),X(2*em_estacao(j)),X(2*de_estacao(j)-1),X(2*de_estacao(j)),X(2*para_estacao(j)-1),X(2*para_estacao(j)));
                    A(j,2*em_estacao(j)-1)=func(1,1);
                    A(j,2*em_estacao(j))=func(2,1);
                    A(j,2*de_estacao(j)-1)=func(3,1);
                    A(j,2*de_estacao(j))=func(4,1);
                    A(j,2*para_estacao(j)-1)=func(5,1);
                    A(j,2*para_estacao(j))=func(6,1);
                    j=j+1;
                end
                % disth
                for j=n_angaz+1:n_angaz+n_disth
                    func=dist_hor_lin(X(2*de_distancia(j-n_angaz)-1),X(2*de_distancia(j-n_angaz)),X(2*para_distancia(j-n_angaz)-1),X(2*para_distancia(j-n_angaz)));
                    A(j,2*de_distancia(j-n_angaz)-1)=func(1,1);
                    A(j,2*de_distancia(j-n_angaz))=func(2,1);
                    A(j,2*para_distancia(j-n_angaz)-1)=func(3,1);
                    A(j,2*para_distancia(j-n_angaz))=func(4,1);
                    j=j+1; 
                end
                % rumo
                for j=n_angaz+n_disth+1:n_angaz+n_disth+n_rumo
                    func=rumo_lin(X(2*de_rumo(j-n_angaz-n_disth)-1),X(2*de_rumo(j-n_angaz-n_disth)),X(2*para_rumo(j-n_angaz-n_disth)-1),X(2*para_rumo(j-n_angaz-n_disth)));
                    A(j,2*de_rumo(j-n_angaz-n_disth)-1)=func(1,1);
                    A(j,2*de_rumo(j-n_angaz-n_disth))=func(2,1);
                    A(j,2*para_rumo(j-n_angaz-n_disth)-1)=func(3,1);
                    A(j,2*para_rumo(j-n_angaz-n_disth))=func(4,1);
                    j=j+1;
                end
                % diraz   
                for j=n_angaz+n_disth+n_rumo+1:total_obs
                    func=direccao_lin(X(2*em_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*em_direccao(j-n_angaz-n_disth-n_rumo)),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)));
                    A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(1,1);
                    A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo))=func(2,1);
                    A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(3,1);
                    A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo))=func(4,1);
                    j=j+1;
                end
                % colunas da matriz A que contem o coeficiente do R0 em cada estacao
                k=1;
                j1=1;
                for j=2:n_diraz
                    if (em_direccao(j) == em_direccao(j-1))
                        A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1; 
                    else
                        A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1;
                        k=k+1;
                    end 
                    j=j+1;
                end
                if (n_diraz>0)
                    A(n_angaz+n_disth+n_rumo+n_diraz,2*n_mp+k)=1;
                end
            end
       
% para eliminar uma coluna, utilizar o operador empty []

            A1=A;

            corta=[c1 c2];
            
            A1(:,corta)=[];

            A1
    
            A2=A(:,[c1 c2]);
            
            A2
     
% calculo da matriz L tal que A1*L=A2 (A1_1 e a inversa generalizada de A1)

            A1_1=pinv(A1);

            L=A1_1*A2;

% separacao do vector X em X1 e X2, X2 sendo o datum

            X1=X;

            X1T=transpose(X1);

            X1T(:,corta)=[];
            
            X1=transpose(X1T)
            
            XT=transpose(X);
            
            X2T=[XT(:,[c1 c2])];
            
            X2=transpose (X2T)
    
% calculo dos angulos utilizando as coordenadas aproximadas

            for j=1:n_angaz
                angaz_c(j)=ang_az_cal(X(2*em_estacao(j)-1),X(2*em_estacao(j)),X(2*de_estacao(j)-1),X(2*de_estacao(j)),X(2*para_estacao(j)-1),X(2*para_estacao(j)));
                j=j+1;
            end
    
            angaz_c=angaz_c';

% calculo das distancias utilizando as coordenadas aproximadas

            for j=1:n_disth
                dist_c(j)=dist_hor_cal(X(2*de_distancia(j)-1),X(2*de_distancia(j)),X(2*para_distancia(j)-1),X(2*para_distancia(j)));
                j=j+1;
            end
            
% calculo dos rumos utilizando as coordenadas aproximadas

            for j=1:n_rumo
                rumo_c=rumo_cal(X(2*de_rumo(j)-1),X(2*de_rumo(j)),X(2*para_rumo(j)-1),X(2*para_rumo(j)));
                j=j+1;
            end 
                
% calculo do vector W utilizando oa angulos observados e os angulos
% calculados    ESTA MAL

            for j=1:n_angaz
                W(j,1)=angaz(j)-angaz_c(j);
                if (W(j,1)<-pi);
                    W(j,1)=W(j,1)+2*pi;
                end        
                if (W(j,1)>pi);
                    W(j,1)=W(j,1)-2*pi;
                end
                j=j+1;
            end
            for j=n_angaz+1:n_angaz+n_disth
                W(j,1)=dist(j-n_angaz)-dist_c(j-n_angaz);
                j=j+1;
            end
            for j=n_angaz+n_disth+1:n_angaz+n_disth+n_rumo
                W(j,1)=rumo(j-n_angaz-n_disth)-rumo_c(j-n_angaz-n_disth);
                j=j+1;
            end
    
            W

% calculo do acrescimo as coordenadas aproximadas

            dX=inv(A1'*Pesos*A1)*A1'*Pesos*W;
        
            dX
           
% calculo das novas coordenadas aproximadas

            if (max(abs(dX)) > conv)
                j=1;
                for i=1:2*n_mp; 
                    if (i~=c1 & i~=c2)
                        X(i,1)=X1(j,1)+dX(j,1);
                        j=j+1;
                    end
                    i=i+1;
                end
            else
                disp ('O processo iterativo convergiu ...')
                break;
            end
            
            iteracao=iteracao+1;
            
            if (iteracao>maxit)
                disp ('O processo iterativo nao convergiu ...')                
                stop
            end  
            
        end

        disp 'Vector dos parametros ajustados:'
        X
        
        QX1=inv((A1'*Pesos*A1));

        ii=1;
        for i=1:2*n_mp;
            if (i~=c1 & i~=c2)
                jj=1;
                for j=1:2*n_mp;
                    if (j~=c1 & j~=c2)
                        QX(i,j)=QX1(ii,jj);
                        jj=jj+1;
                    end
                    j=j+1;
                end
                ii=ii+1;
            end
            i=i+1;
        end
        
%        disp 'Matriz cofactor das coordenadas ajustadas:'
        QX;
        
        tracoQX=0;
        for i=1:2*n_mp;
            tracoQX=tracoQX+QX(i,i);
            i=i+1;
        end
        
        disp 'Traco da matriz cofactor das coordenadas ajustadas:'
        tracoQX
        
        disp 'Vector dos residuos das observacoes:'        
        v=A1*dX+W
        
        Qv=Q-A1*inv(A1'*Pesos*A1)*A1';
        
        l_adj=l+v;
        for i=1:n_angaz
            l_adj=l_adj*180/pi;
        end
        
        l_adj;
        
        Ql_adj=Q-Qv;
        
        sigma_posteriori_2=v'*Pesos*v/gl;
        
        disp 'Teste da razao de variancias:'
        disp '                             '
 
        qui=icdf('chi2',1-nsign,gl);
        
        if (sigma_posteriori_2*gl/sigma_priori_2 < qui)
            disp 'Aceitar o ajustamento'
        else
            disp 'Rejeitar o ajustamento'
        end
        
    end
       
end

% *************************************************************************
% ******************* OPCAO = 2 *******************************************
% *************************************************************************

if (opcao1 == 2)

    clc
    disp '**************************************************************'
    disp '*                 Datum com constrangimentos                 *'
    disp '**************************************************************'    
    disp '                                                              '

    if (n_diraz>0)
        Rt=zeros(nulidade,2*n_mp+ne);  
    else
        Rt=zeros(nulidade,2*n_mp);
    end    
    c=zeros(nulidade,1);
    opcao2=0;
    
    disp ('Tipos de constrangimento numa rede 2D ');
    disp '                                                          '
    disp ('  1. Fixar as coordenadas (X,Y) de 2 pontos da rede (d-4)');
    disp ('  2. Fixar as coordenadas (X,Y) de 1 ponto da rede (d-2)');
    disp ('  3. Fixar a coordenada X de 1 ponto da rede (d-1)');
    disp ('  4. Fixar a coordenada Y de 1 ponto da rede (d-1)');
    disp ('  5. Fixar a distancia horizontal entre 2 pontos de rede (d-1)');
    disp ('  6. Fixar o rumo da direccao entre 2 pontos de rede (d-1)');
    disp ('  7. Fixar a diferenca de coordenadas (X,Y) entre 2 pontos de rede (d-2)');
    disp '                                                          '

% calculo da matriz dos constrangimentos em funcao da opcao escolhida    

% *************************************************************************

    if (n_disth == 0 & n_rumo == 0) %nulidade == 4: direccoes ou angulos azimutais
        disp(['Opcoes possiveis para o caso de o deficit de caracteristica da matriz N ser igual a d = ',num2str(nulidade)]); 
        disp ('                                                    ');
        disp ('                         1')
        disp ('               256 (2+5+6)')
        disp ('                                                    ');   
        opcao2=input('     Seleccione a opcao: ');
    end
    
    clc
                                                     
    if (n_disth == 0 & n_rumo == 0 & opcao2 == 1) % 2 pontos fixos
        
        disp ('                                                    ');   
        for j=1:n_mp;
            if (j<10)
                str=['        ',num2str(j),' - ',char(ponto(j))];
            else
                str=['       ',num2str(j),' - ',char(ponto(j))];        
            end
            disp (str)
            j=j+1;
        end        
        disp ('                                                     ');
        datum1=input ('  Indique o primeiro ponto seleccionado para datum:  ');
        
        disp ('                                                    ');   
        for j=1:n_mp;
            if (j ~=datum1)
                if (j<10)
                    str=['        ',num2str(j),' - ',char(ponto(j))];
                else
                    str=['       ',num2str(j),' - ',char(ponto(j))];        
                end
                disp (str)
                j=j+1;
            end
        end      
        disp ('                                                    ');   
        datum2=input ('  Indique  o segundo ponto seleccionado para datum:  ');
        
        Rt(1,2*datum1-1)=1;
        Rt(2,2*datum1)=1;
        Rt(3,2*datum2-1)=1;
        Rt(4,2*datum2)=1;
    end

    if (n_disth == 0 & n_rumo == 0 & opcao2 == 256) % 1 ponto fixo + 1 distancia + 1 rumo
        disp ('                                                    ');   
        for j=1:n_mp;
            if (j<10)
                str=['        ',num2str(j),' - ',char(ponto(j))];
            else
                str=['       ',num2str(j),' - ',char(ponto(j))];        
            end
            disp (str)
            j=j+1;
        end        
        disp ('                                                    ');
        datum1=input ('  Indique o ponto seleccionado para datum:  ');
        disp ('                                                    ');
        disp ('  Fixe a distancia entre: ');
        ponto1=input ('    Ponto inicial: ');
        ponto2=input ('      Ponto final: ');
        disp ('                                                    ');
        disp ('  Fixe o rumo entre: ');
        ponto3=input ('    Ponto inicial: ');
        ponto4=input ('      Ponto final: ');
        Rt(1,2*datum1-1)=1;
        Rt(2,2*datum1)=1;
        Rt(3,2*ponto1-1)=-(X(2*ponto2-1)-X(2*ponto1-1));
        Rt(3,2*ponto1)=-(X(2*ponto2)-X(2*ponto1));
        Rt(3,2*ponto2-1)=(X(2*ponto2-1)-X(2*ponto1-1));
        Rt(3,2*ponto2)=(X(2*ponto2)-X(2*ponto1));
        Rt(4,2*ponto3-1)=-(X(2*ponto4)-X(2*ponto3));
        Rt(4,2*ponto3)=(X(2*ponto4-1)-X(2*ponto3-1));
        Rt(4,2*ponto4-1)=(X(2*ponto4)-X(2*ponto3));
        Rt(4,2*ponto4)=-(X(2*ponto4-1)-X(2*ponto3-1));
    end    

% *************************************************************************

    if (n_disth > 0 & n_rumo == 0) %nulidade == 3
        disp(['Opcoes possiveis para o caso de o deficit de caracteristica da matriz N ser igual a d = ',num2str(nulidade),' (com obs. de distancias)']); 
        disp ('                                                    ');
        disp ('                 23 (2+3)')
        disp ('                 24 (2+4)')
        disp ('                 26 (2+6)')    
        disp ('                                                    ');   
        opcao2=input('     Seleccione a opcao: ');
    end

    if (n_disth > 0 & n_rumo == 0 & opcao2==23) % 1 ponto fixo + X de outro ponto
        disp ('                                                    ');   
        for j=1:n_mp;
            if (j<10)
                str=['        ',num2str(j),' - ',char(ponto(j))];
            else
                str=['       ',num2str(j),' - ',char(ponto(j))];        
            end
            disp (str)
            j=j+1;
        end                
        disp ('                                                    ');           
        datum1=input ('  Indique o primeiro ponto seleccionado para datum:  ');

        disp '                                                          '
        for j=1:n_mp;
            if (j ~= datum1)
                if (j<10)
                    str=['        ',num2str(j),' - X',char(ponto(j))];
                else
                    str=['       ',num2str(j),' - X',char(ponto(j))];        
                end
                disp (str)
                j=j+1;
            end
        end                
        disp ('                                                         ');       
        datum2=input ('  Indique o segundo ponto seleccionado para datum (coordenada X):  ');     
        
        Rt(1,2*datum1-1)=1;
        Rt(2,2*datum1)=1;
        Rt(3,2*datum2-1)=1;   
    end

    if (n_disth > 0 & n_rumo == 0 & opcao2==24) % 1 ponto fixo + Y de outro ponto
        disp ('                                                    ');   
        for j=1:n_mp;
            if (j<10)
                str=['        ',num2str(j),' - ',char(ponto(j))];
            else
                str=['       ',num2str(j),' - ',char(ponto(j))];        
            end
            disp (str)
            j=j+1;
        end                
        disp ('                                                    ');           
        datum1=input ('  Indique o primeiro ponto seleccionado para datum:  ');

        disp '                                                          '
        for j=1:n_mp;
            if (j ~= datum1)
                if (j<10)
                    str=['        ',num2str(j),' - Y',char(ponto(j))];
                else
                    str=['       ',num2str(j),' - Y',char(ponto(j))];        
                end
                disp (str)
                j=j+1;
            end
        end                
        disp ('                                                         ');         
        datum2=input ('  Indique o segundo ponto seleccionado para datum (coordenada Y):  ');
        
        Rt(1,2*datum1-1)=1;
        Rt(2,2*datum1)=1;
        Rt(3,2*datum2)=1;
    end

    if (n_disth > 0 & n_rumo == 0 & opcao2==26) % 1 ponto fixo + 1 rumo
        disp ('                                                    ');   
        for j=1:n_mp;
            if (j<10)
                str=['        ',num2str(j),' - ',char(ponto(j))];
            else
                str=['       ',num2str(j),' - ',char(ponto(j))];        
            end
            disp (str)
            j=j+1;
        end        
        disp ('                                                    ');                      
        datum1=input ('  Indique o ponto seleccionado para datum:  ');
        
        disp ('                                                    ');                      
        disp ('  Fixe o rumo entre: ');
        ponto3=input ('    Ponto inicial: ');
        ponto4=input ('      Ponto final: ');
        
        Rt(1,2*datum1-1)=1;
        Rt(2,2*datum1)=1;
        Rt(3,2*ponto3-1)=-(X(2*ponto4)-X(2*ponto3));
        Rt(3,2*ponto3)=(X(2*ponto4-1)-X(2*ponto3-1));
        Rt(3,2*ponto4-1)=(X(2*ponto4)-X(2*ponto3));
        Rt(3,2*ponto4)=-(X(2*ponto4-1)-X(2*ponto3-1));
    end

    if (n_rumo > 0 & n_disth == 0) %nulidade == 3
        disp(['Opcoes possiveis para o caso de o deficit de caracteristica da matriz N ser igual a d = ',num2str(nulidade),' (com obs. de rumo)']); 
        disp ('                                                    ');
        disp ('                      23')
        disp ('                      24')
        disp ('                      25')   
        disp ('                                                    ');   
        opcao2=input('     Seleccione a opcao: ');
    end

     if (n_rumo > 0 & n_disth == 0 & opcao2==23) % 1 ponto fixo + X de outro ponto 
        disp ('                                                    ');   
        for j=1:n_mp;
            if (j<10)
                str=['        ',num2str(j),' - ',char(ponto(j))];
            else
                str=['       ',num2str(j),' - ',char(ponto(j))];        
            end
            disp (str)
            j=j+1;
        end                
        disp ('                                                    ');           
        datum1=input ('  Indique o primeiro ponto seleccionado para datum:  ');
        
        disp '                                                          '
        for j=1:n_mp;
            if (j ~= datum1)
                if (j<10)
                    str=['        ',num2str(j),' - X',char(ponto(j))];
                else
                    str=['       ',num2str(j),' - X',char(ponto(j))];        
                end
                disp (str)
                j=j+1;
            end
        end                
        disp ('                                                         ');       
        datum2=input ('  Indique o segundo ponto seleccionado para datum (coordenada X):  ');                
        
        Rt(1,2*datum1-1)=1;
        Rt(2,2*datum1)=1;
        Rt(3,2*datum2-1)=1;
     end

    
     if (n_disth == 0 & n_rumo > 0 & opcao2==24) % 1 ponto fixo + Y de outro ponto
        disp ('                                                    ');   
        for j=1:n_mp;
            if (j<10)
                str=['        ',num2str(j),' - ',char(ponto(j))];
            else
                str=['       ',num2str(j),' - ',char(ponto(j))];        
            end
            disp (str)
            j=j+1;
        end                
        disp ('                                                    ');           
        datum1=input ('  Indique o primeiro ponto seleccionado para datum:  ');

        disp '                                                          '
        for j=1:n_mp;
            if (j ~= datum1)
                if (j<10)
                    str=['        ',num2str(j),' - Y',char(ponto(j))];
                else
                    str=['       ',num2str(j),' - Y',char(ponto(j))];        
                end
                disp (str)
                j=j+1;
            end
        end                
        disp ('                                                         ');   
        datum2=input ('  Indique o segundo ponto seleccionado para datum (coordenada X):  ');   
        
        Rt(1,2*datum1-1)=1;
        Rt(2,2*datum1)=1;
        Rt(3,2*datum2-1)=1;   
    end
   
    if (n_rumo > 0 & n_disth == 0 & opcao2==25) % 1 ponto fixo + 1 distancia
        disp ('                                                    ');   
        for j=1:n_mp;
            if (j<10)
                str=['        ',num2str(j),' - ',char(ponto(j))];
            else
                str=['       ',num2str(j),' - ',char(ponto(j))];        
            end
            disp (str)
            j=j+1;
        end   
        disp ('                                                    ');
        datum1=input ('  Indique o ponto seleccionado para datum:  ');
        disp ('  Fixe a distancia entre: ');
        ponto1=input ('    Ponto inicial: ');
        ponto2=input ('      Ponto final: ');
        Rt(1,2*datum1-1)=1;
        Rt(2,2*datum1)=1;
        Rt(3,2*ponto1-1)=-(X(2*ponto2-1)-X(2*ponto1-1));
        Rt(3,2*ponto1)=-(X(2*ponto2)-X(2*ponto1));
        Rt(3,2*ponto2-1)=(X(2*ponto2-1)-X(2*ponto1-1));
        Rt(3,2*ponto2)=(X(2*ponto2)-X(2*ponto1));
    end

% *************************************************************************

    if (n_disth > 0 & n_rumo > 0) % nulidade == 2
        disp(['Opcoes possiveis para o caso de o deficit de caracteristica da matriz N ser igual a d = ',num2str(nulidade)]); 
        disp ('                                                    ');
        disp ('                       2')
        opcao2=input('     Seleccione a opcao: ');
    end

    if (n_disth > 0 & n_rumo > 0 & opcao2==2) % 1 ponto fixo
        disp ('                                                    ');   
        for j=1:n_mp;
            if (j<10)
                str=['        ',num2str(j),' - ',char(ponto(j))];
            else
                str=['       ',num2str(j),' - ',char(ponto(j))];        
            end
            disp (str)
            j=j+1;
        end    
        disp ('                                                    ');
        datum1=input ('  Indique o ponto seleccionado para datum:  ');
        Rt(1,2*datum1-1)=1;
        Rt(2,2*datum1)=1;
    end

    Rt;

    R=Rt';

% inicio do processo iterativo

    for iteracao=1:maxit
    
        iteracao

% actualizacao da matriz A-----------

        if (iteracao>1)
            % angaz
            for j=1:n_angaz;
                func=ang_az_lin(X(2*em_estacao(j)-1),X(2*em_estacao(j)),X(2*de_estacao(j)-1),X(2*de_estacao(j)),X(2*para_estacao(j)-1),X(2*para_estacao(j)));
                A(j,2*em_estacao(j)-1)=func(1,1);
                A(j,2*em_estacao(j))=func(2,1);
                A(j,2*de_estacao(j)-1)=func(3,1);
                A(j,2*de_estacao(j))=func(4,1);
                A(j,2*para_estacao(j)-1)=func(5,1);
                A(j,2*para_estacao(j))=func(6,1);
                j=j+1;
            end
            % disth
            for j=n_angaz+1:n_angaz+n_disth
                func=dist_hor_lin(X(2*de_distancia(j-n_angaz)-1),X(2*de_distancia(j-n_angaz)),X(2*para_distancia(j-n_angaz)-1),X(2*para_distancia(j-n_angaz)));
                A(j,2*de_distancia(j-n_angaz)-1)=func(1,1);
                A(j,2*de_distancia(j-n_angaz))=func(2,1);
                A(j,2*para_distancia(j-n_angaz)-1)=func(3,1);
                A(j,2*para_distancia(j-n_angaz))=func(4,1);
                j=j+1; 
            end
            % rumo
            for j=n_angaz+n_disth+1:n_angaz+n_disth+n_rumo
                func=rumo_lin(X(2*de_rumo(j-n_angaz-n_disth)-1),X(2*de_rumo(j-n_angaz-n_disth)),X(2*para_rumo(j-n_angaz-n_disth)-1),X(2*para_rumo(j-n_angaz-n_disth)));
                A(j,2*de_rumo(j-n_angaz-n_disth)-1)=func(1,1);
                A(j,2*de_rumo(j-n_angaz-n_disth))=func(2,1);
                A(j,2*para_rumo(j-n_angaz-n_disth)-1)=func(3,1);
                A(j,2*para_rumo(j-n_angaz-n_disth))=func(4,1);
                j=j+1;
            end
            % diraz   
            for j=n_angaz+n_disth+n_rumo+1:total_obs
                func=direccao_lin(X(2*em_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*em_direccao(j-n_angaz-n_disth-n_rumo)),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)));
                A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(1,1);
                A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo))=func(2,1);
                A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(3,1);
                A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo))=func(4,1);
                j=j+1;
            end
            % colunas da matriz A que contem o coeficiente do R0 em cada estacao
            k=1;
            j1=1;
            for j=2:n_diraz
                if (em_direccao(j) == em_direccao(j-1))
                    A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1; 
                else
                    A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1;
                    k=k+1;
                end 
                j=j+1;
            end
            if (n_diraz>0)
                A(n_angaz+n_disth+n_rumo+n_diraz,2*n_mp+k)=1;
            end
        end
   
        A;
        
% actualizacao da matriz R
        
        if (iteracao>1)
            if (n_disth == 0 & n_rumo == 0) %nulidade == 4
                if (opcao2 == 1) % 2 pontos fixos
                    Rt(1,2*datum1-1)=1;
                    Rt(2,2*datum1)=1;
                    Rt(3,2*datum2-1)=1;
                    Rt(4,2*datum2)=1;
                end
                if (opcao2 == 256) % 1 ponto fixo + 1 distancia + 1 rumo
                    Rt(1,2*datum1-1)=1;
                    Rt(2,2*datum1)=1;
                    Rt(3,2*ponto1-1)=-(X(2*ponto2-1)-X(2*ponto1-1));
                    Rt(3,2*ponto1)=-(X(2*ponto2)-X(2*ponto1));
                    Rt(3,2*ponto2-1)=(X(2*ponto2-1)-X(2*ponto1-1));
                    Rt(3,2*ponto2)=(X(2*ponto2)-X(2*ponto1));
                    Rt(4,2*ponto3-1)=-(X(2*ponto4)-X(2*ponto3));
                    Rt(4,2*ponto3)=(X(2*ponto4-1)-X(2*ponto3-1));
                    Rt(4,2*ponto4-1)=(X(2*ponto4)-X(2*ponto3));
                    Rt(4,2*ponto4)=-(X(2*ponto4-1)-X(2*ponto3-1));
                end                    
            end                                                    
            if (n_disth > 0 & n_rumo == 0) %nulidade == 3
                if (n_disth > 0 & n_rumo == 0 & opcao2==23) % 1 ponto fixo + X de outro ponto                                                           ');
                    Rt(1,2*datum1-1)=1;
                    Rt(2,2*datum1)=1;
                    Rt(3,2*datum2-1)=1;   
                end        
                if (n_disth > 0 & n_rumo == 0 & opcao2==24) % 1 ponto fixo + Y de outro ponto
                    Rt(1,2*datum1-1)=1;
                    Rt(2,2*datum1)=1;
                    Rt(3,2*datum2)=1;
                end           
                if (opcao2==26) % 1 ponto fixo + 1 rumo
                    Rt(1,2*datum1-1)=1;
                    Rt(2,2*datum1)=1;
                    Rt(3,2*ponto3-1)=-(X(2*ponto4)-X(2*ponto3));
                    Rt(3,2*ponto3)=(X(2*ponto4-1)-X(2*ponto3-1));
                    Rt(3,2*ponto4-1)=(X(2*ponto4)-X(2*ponto3));
                    Rt(3,2*ponto4)=-(X(2*ponto4-1)-X(2*ponto3-1));
                end              
            end
            if (n_rumo > 0 & n_disth == 0) %nulidade == 3
                if (n_rumo > 0 & n_disth == 0 & opcao2==23) % 1 ponto fixo + X de outro ponto
                    Rt(1,2*datum1-1)=1;
                    Rt(2,2*datum1)=1;
                    Rt(3,2*datum2-1)=1;   
                end
                if (n_rumo > 0 & n_disth == 0 & opcao2==24) % 1 ponto fixo + Y de outro ponto
                    Rt(1,2*datum1-1)=1;
                    Rt(2,2*datum1)=1;
                    Rt(3,2*datum2)=1;
                end     
                if (opcao2==25) % 1 ponto fixo + 1 distancia
                    Rt(1,2*datum1-1)=1;
                    Rt(2,2*datum1)=1;
                    Rt(3,2*ponto1-1)=-(X(2*ponto2-1)-X(2*ponto1-1));
                    Rt(3,2*ponto1)=-(X(2*ponto2)-X(2*ponto1));
                    Rt(3,2*ponto2-1)=(X(2*ponto2-1)-X(2*ponto1-1));
                    Rt(3,2*ponto2)=(X(2*ponto2)-X(2*ponto1));
                end                
            end
            if (n_disth > 0 & n_rumo > 0) % nulidade == 2
                if (opcao2==2) % 1 ponto fixo
                    Rt(1,2*datum1-1)=1;
                    Rt(2,2*datum1)=1;
                end          
            end
        end
         
% calculo do vector fecho W:

% calculo das distancias utilizando as coordenadas aproximadas

        for j=1:n_disth
            dist_c(j)=dist_hor_cal(X(2*de_distancia(j)-1),X(2*de_distancia(j)),X(2*para_distancia(j)-1),X(2*para_distancia(j)));
            j=j+1;
        end
        
        if (n_disth>0)
            dist_c=dist_c';
        end
              
% transformacao de direccoes azimutais em rumos atraves do calculo do R0 em
% cada estacao; para o efeito, e necessario contar o numero de estacoes
% ocupadas e em cada estacao o numero de pontos visados; em  cada estacao
% calcula-se o rumo ponto estacao-ponto visado, tendo-se entao R0=REV-LEV

% calculo do R0 para cada leitura azimutal
    
        for j=1:n_diraz
            R0(j)=rumo(M(em_direccao(j)),P(em_direccao(j)),M(para_direccao(j)),P(para_direccao(j)))-diraz(j);
            if (R0(j)<0)
                R0(j)=R0(j)+2*pi;
            end
        end

% calculo do R0 medio e respectivo desvio padrao em cada estacao 

        sum_rumo=R0(1);
        sum_rumo2=sum_rumo^2;
        num_pontarias=1;
        j1=1;
        for j=2:n_diraz
            if (em_direccao(j) == em_direccao(j-1))
                sum_rumo=sum_rumo+R0(j);
                sum_rumo2=sum_rumo2+R0(j)^2;
                num_pontarias=num_pontarias+1;
            else
                j2=j-1;
                rumo_medio=sum_rumo/num_pontarias;        
                std_rumo_medio=sqrt((sum_rumo2-sum_rumo^2/(num_pontarias))/(num_pontarias-1));
                for k=j1:j2 
                    R0_medio(k)=rumo_medio;
                    sigma_R0_medio(k)=std_rumo_medio;
                    k=k+1;
                end
                sum_rumo=R0(j);
                sum_rumo2=sum_rumo^2;
                num_pontarias=1;
                j1=j;
            end 
            j=j+1;
        end

        j2=j-1;
        rumo_medio=sum_rumo/num_pontarias;        
        std_rumo_medio=sqrt((sum_rumo2-sum_rumo^2/(num_pontarias))/(num_pontarias-1));  
        for k=j1:j2 
            R0_medio(k)=rumo_medio;
            sigma_R0_medio(k)=std_rumo_medio;
            k=k+1;
        end

        R0=R0';
        R0_medio=R0_medio';
        sigma_R0_medio=sigma_R0_medio';

% calculo das direccoes azimutais "calculadas", dadas como DIRAZ=R-R0

        j1=1;
        for j=2:n_diraz
            if (em_direccao(j) == em_direccao(j-1))
                diraz_c(j-1)=rumo_cal(X(2*em_direccao(j-1)-1),X(2*em_direccao(j-1)),X(2*para_direccao(j-1)-1),X(2*para_direccao(j-1)))-R0_medio(j1);
                if (diraz_c(j-1)<0)
                    diraz_c(j-1)=diraz_c(j-1)+2*pi;
                end
            else
               diraz_c(j-1)=rumo_cal(X(2*em_direccao(j-1)-1),X(2*em_direccao(j-1)),X(2*para_direccao(j-1)-1),X(2*para_direccao(j-1)))-R0_medio(j1); 
               if (diraz_c(j-1)<0)
                   diraz_c(j-1)=diraz_c(j-1)+2*pi;
               end                   
               j1=j;
            end 
            j=j+1;
        end
        diraz_c(n_diraz)=rumo_cal(X(2*em_direccao(n_diraz)-1),X(2*em_direccao(n_diraz)),X(2*para_direccao(n_diraz)-1),X(2*para_direccao(n_diraz)))-R0_medio(j1); 
        if (diraz_c(n_diraz)<0)
            diraz_c(n_diraz)=diraz_c(n_diraz)+2*pi;
        end     
                   
        diraz_c=diraz_c';
        
% calculo dos rumos utilizando as coordenadas aproximadas

        for j=1:n_rumo
            rumo_c=rumo_cal(X(2*de_rumo(j)-1),X(2*de_rumo(j)),X(2*para_rumo(j)-1),X(2*para_rumo(j)));
            j=j+1;
        end         
                      
% calculo do vector W (vector fecho)

        for j=1:n_angaz
            W(j,1)=angaz_c(j)-angaz(j);
            if (W(j,1)<-pi);
                W(j,1)=W(j,1)+2*pi;
            end        
            if (W(j,1)>pi);
                W(j,1)=W(j,1)-2*pi;
            end
            j=j+1;
        end           
        for j=n_angaz+1:n_angaz+n_disth
            W(j,1)=dist_c(j-n_angaz)-disth(j-n_angaz);
            j=j+1;
        end
        for j=n_angaz+n_disth+1:n_angaz+n_disth+n_rumo
            W(j,1)=rumo(j-n_angaz-n_disth)-rumo_c(j-n_angaz-n_disth);
            j=j+1;
        end        
        for j=1:n_diraz
            W(n_disth+j,1)=diraz_c(j)-diraz(j);
            if (W(n_disth+j,1)<-pi);
                W(n_disth+j,1)=W(n_disth+j,1)+2*pi;
            end        
            if (W(n_disth+j,1)>pi);
                W(n_disth+j,1)=W(n_disth+j,1)-2*pi;
            end
            j=j+1;
        end
    
        W;
        
% calculo do acrescimo para as coordenadas aproximadas

        dX=-inv(N+R*Rt)*N*inv(N+R*Rt)*A'*Pesos*W + inv(N+R*Rt)*R*c;

% calculo das novas coordenadas aproximadas

        if (max(abs(dX))>conv)
%%            for i=1:2*n_mp+ne;
            for i=1:2*n_mp;                
                X(i,1)=X(i,1)+dX(i,1);
                i=i+1;
            end
        else
            disp ('O processo iterativo convergiu.')
            disp '                                  '
            break;
        end

        if (iteracao > maxit)
            disp ('O processo iterativo nao convergiu ...')
            stop
        end
        
        iteracao=iteracao+1;
        
    end

    disp 'Vector das coordenadas ajustadas:´'
    X

    QX=inv(N+R*Rt)*N*inv(N+R*Rt)

    tracoQX=0;
    for i=1:2*n_mp
        tracoQX=tracoQX+QX(i,i);
    end

    disp 'Traco da matriz cofactor das coordenadas ajustadas'
    tracoQX

    disp 'Vector dos residuos das observacoes:'
    v=A*dX+W
    
    Qv=Ql-A*inv(A'*Pesos*A)*A';

    l_adj=l+v;
    for i=1:n_angaz
        l_adj(i,1)=l_adj(i,1)*200/pi;
        i=i+1;
    end
    for i=n_angaz+n_rumo+n_disth+1:total_obs
        l_adj(i,1)=l_adj(i,1)*200/pi;
        i=i+1;
    end
        
    l_adj;

    Ql_adj=Ql-Qv;

    sigma_posteriori_2=v'*Pesos*v/gl;
        
    disp 'Teste da razao de variancias:'
    disp '                             '
 
    qui=icdf('chi2',1-nsign,gl);
        
    if (sigma_posteriori_2*gl/sigma_priori_2 < qui)
        disp 'Aceitar o ajustamento'
    else
        disp 'Rejeitar o ajustamento'
    end

% elipses de erro e/ou elipses de confianca    
    
end

% *************************************************************************
% ******************* OPCAO = 3 *******************************************
% *************************************************************************

if (opcao1 == 3)

    clc
    disp '**************************************************************'
    disp '*                 Datum com traco minimo                     *'
    disp '**************************************************************'    
    disp '                                                              '
    
    
    for iteracao=1:maxit
    
        iteracao

% actualizacao da matriz A

        if (iteracao>1)
            % angaz
            for j=1:n_angaz;
                func=ang_az_lin(X(2*em_estacao(j)-1),X(2*em_estacao(j)),X(2*de_estacao(j)-1),X(2*de_estacao(j)),X(2*para_estacao(j)-1),X(2*para_estacao(j)));
                A(j,2*em_estacao(j)-1)=func(1,1);
                A(j,2*em_estacao(j))=func(2,1);
                A(j,2*de_estacao(j)-1)=func(3,1);
                A(j,2*de_estacao(j))=func(4,1);
                A(j,2*para_estacao(j)-1)=func(5,1);
                A(j,2*para_estacao(j))=func(6,1);
                j=j+1;
            end
            % disth
            for j=n_angaz+1:n_angaz+n_disth
                func=dist_hor_lin(X(2*de_distancia(j-n_angaz)-1),X(2*de_distancia(j-n_angaz)),X(2*para_distancia(j-n_angaz)-1),X(2*para_distancia(j-n_angaz)));
                A(j,2*de_distancia(j-n_angaz)-1)=func(1,1);
                A(j,2*de_distancia(j-n_angaz))=func(2,1);
                A(j,2*para_distancia(j-n_angaz)-1)=func(3,1);
                A(j,2*para_distancia(j-n_angaz))=func(4,1);
                j=j+1; 
            end
            % rumo
            for j=n_angaz+n_disth+1:n_angaz+n_disth+n_rumo
                func=rumo_lin(X(2*de_rumo(j-n_angaz-n_disth)-1),X(2*de_rumo(j-n_angaz-n_disth)),X(2*para_rumo(j-n_angaz-n_disth)-1),X(2*para_rumo(j-n_angaz-n_disth)));
                A(j,2*de_rumo(j-n_angaz-n_disth)-1)=func(1,1);
                A(j,2*de_rumo(j-n_angaz-n_disth))=func(2,1);
                A(j,2*para_rumo(j-n_angaz-n_disth)-1)=func(3,1);
                A(j,2*para_rumo(j-n_angaz-n_disth))=func(4,1);
                j=j+1;
            end
            % diraz   
            for j=n_angaz+n_disth+n_rumo+1:total_obs
                func=direccao_lin(X(2*em_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*em_direccao(j-n_angaz-n_disth-n_rumo)),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)));
                A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(1,1);
                A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo))=func(2,1);
                A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(3,1);
                A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo))=func(4,1);
                j=j+1;
            end
            % colunas da matriz A que contem o coeficiente do R0 em cada estacao
            k=1;
            j1=1;
            for j=2:n_diraz
                if (em_direccao(j) == em_direccao(j-1))
                    A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1; 
                else
                    A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1;
                    k=k+1;
                end 
                j=j+1;
            end
            if (n_diraz>0)
                A(n_angaz+n_disth+n_rumo+n_diraz,2*n_mp+k)=1;
            end
        end
                
% calculo do vector fecho W:

% calculo das distancias utilizando as coordenadas aproximadas

        for j=1:n_disth
            dist_c(j)=dist_hor_cal(X(2*de_distancia(j)-1),X(2*de_distancia(j)),X(2*para_distancia(j)-1),X(2*para_distancia(j)));
            j=j+1;
        end
        
        if (n_disth>0)
            dist_c=dist_c';
        end
              
% transformacao de direccoes azimutais em rumos atraves do calculo do R0 em
% cada estacao; para o efeito, e necessario contar o numero de estacoes
% ocupadas e em cada estacao o numero de pontos visados; em  cada estacao
% calcula-se o rumo ponto estacao-ponto visado, tendo-se entao R0=REV-LEV

% calculo do R0 para cada leitura azimutal
    
        for j=1:n_diraz
            R0(j)=rumo(M(em_direccao(j)),P(em_direccao(j)),M(para_direccao(j)),P(para_direccao(j)))-diraz(j);
            if (R0(j)<0)
                R0(j)=R0(j)+2*pi;
            end
        end

% calculo do R0 medio e respectivo desvio padrao em cada estacao 

        sum_rumo=R0(1);
        sum_rumo2=sum_rumo^2;
        num_pontarias=1;
        j1=1;
        for j=2:n_diraz
            if (em_direccao(j) == em_direccao(j-1))
                sum_rumo=sum_rumo+R0(j);
                sum_rumo2=sum_rumo2+R0(j)^2;
                num_pontarias=num_pontarias+1;
            else
                j2=j-1;
                rumo_medio=sum_rumo/num_pontarias;        
                std_rumo_medio=sqrt((sum_rumo2-sum_rumo^2/(num_pontarias))/(num_pontarias-1));
                for k=j1:j2 
                    R0_medio(k)=rumo_medio;
                    sigma_R0_medio(k)=std_rumo_medio;
                    k=k+1;
                end
                sum_rumo=R0(j);
                sum_rumo2=sum_rumo^2;
                num_pontarias=1;
                j1=j;
            end 
            j=j+1;
        end

        j2=j-1;
        rumo_medio=sum_rumo/num_pontarias;        
        std_rumo_medio=sqrt((sum_rumo2-sum_rumo^2/(num_pontarias))/(num_pontarias-1));  
        for k=j1:j2 
            R0_medio(k)=rumo_medio;
            sigma_R0_medio(k)=std_rumo_medio;
            k=k+1;
        end

        R0=R0';
        R0_medio=R0_medio';
        sigma_R0_medio=sigma_R0_medio';

% calculo das direccoes azimutais "calculadas", dadas como DIRAZ=R-R0

        j1=1;
        for j=2:n_diraz
            if (em_direccao(j) == em_direccao(j-1))
                diraz_c(j-1)=rumo_cal(X(2*em_direccao(j-1)-1),X(2*em_direccao(j-1)),X(2*para_direccao(j-1)-1),X(2*para_direccao(j-1)))-R0_medio(j1);
                if (diraz_c(j-1)<0)
                    diraz_c(j-1)=diraz_c(j-1)+2*pi;
                end
            else
               diraz_c(j-1)=rumo_cal(X(2*em_direccao(j-1)-1),X(2*em_direccao(j-1)),X(2*para_direccao(j-1)-1),X(2*para_direccao(j-1)))-R0_medio(j1); 
               if (diraz_c(j-1)<0)
                   diraz_c(j-1)=diraz_c(j-1)+2*pi;
               end                   
               j1=j;
            end 
            j=j+1;
        end
        diraz_c(n_diraz)=rumo_cal(X(2*em_direccao(n_diraz)-1),X(2*em_direccao(n_diraz)),X(2*para_direccao(n_diraz)-1),X(2*para_direccao(n_diraz)))-R0_medio(j1); 
        if (diraz_c(n_diraz)<0)
            diraz_c(n_diraz)=diraz_c(n_diraz)+2*pi;
        end     
                   
        diraz_c=diraz_c';
        
% calculo dos rumos utilizando as coordenadas aproximadas

        for j=1:n_rumo
            rumo_c=rumo_cal(X(2*de_rumo(j)-1),X(2*de_rumo(j)),X(2*para_rumo(j)-1),X(2*para_rumo(j)));
            j=j+1;
        end         
                      
% calculo do vector W (vector fecho)

        for j=1:n_angaz
            W(j,1)=angaz_c(j)-angaz(j);
            if (W(j,1)<-pi);
                W(j,1)=W(j,1)+2*pi;
            end        
            if (W(j,1)>pi);
                W(j,1)=W(j,1)-2*pi;
            end
            j=j+1;
        end           
        for j=n_angaz+1:n_angaz+n_disth
            W(j,1)=dist_c(j-n_angaz)-disth(j-n_angaz);
            j=j+1;
        end
        for j=n_angaz+n_disth+1:n_angaz+n_disth+n_rumo
            W(j,1)=rumo(j-n_angaz-n_disth)-rumo_c(j-n_angaz-n_disth);
            j=j+1;
        end        
        for j=1:n_diraz
            W(n_disth+j,1)=diraz_c(j)-diraz(j);
            if (W(n_disth+j,1)<-pi);
                W(n_disth+j,1)=W(n_disth+j,1)+2*pi;
            end        
            if (W(n_disth+j,1)>pi);
                W(n_disth+j,1)=W(n_disth+j,1)-2*pi;
            end
            j=j+1;
        end
    
        W

% calculo do acrescimo as coordenadas aproximadas

% [V,D] = eig(A) produces matrices of eigenvalues (D) and eigenvectors (V) of matrix A, so that A*V = V*D. 
% Matrix D is the canonical form of A - a diagonal matrix with A's eigenvalues on the main diagonal. 
% Matrix V is the modal matrix - its columns are the eigenvectors of A.

        [linhasN,colunasN]=size(N);

        [Vec,Val]=eig(N,'nobalance');
    
        Val;

        Vec;

% matriz dos constrangimentos    
    
        for i=1:linhasN,
            for j=1:linhasN,
                R(i,j)=Vec(i,j);
            end
        end

        R;

% calculo do vector dos parametros

        B=inv(N+R*R')*N*inv(N+R*R');

        dx=-B*A'*Pesos*W
           
% calculo das novas coordenadas aproximadas

        if (max(abs(dx)) > conv)
            for i=1:2*n_mp; 
                X(i,1)=X(i,1)+dx(i,1);
                i=i+1;
            end
        else
            disp ('O processo iterativo convergiu ...')
            break;
        end
        
        if (iteracao > maxit)
            disp ('O processo iterativo nao convergiu ...')
            stop
        end
                    
        iteracao=iteracao+1;
            
    end    
   
    disp 'Vector das coordenadas ajustadas:'    
    X
    
    QX=inv(N+R*R')*N*inv(N+R*R');

    tracoQX=0;
    
    for i=1:2*n_mp
        tracoQX=tracoQX+QX(i,i);
    end
    
    disp 'Traco da matriz cofactor dos parametros ajustados:'
    tracoQX

    disp 'Vector dos residuos ajustados:'
    v=A*dx+W

    Qv=Ql-A*inv(A'*Pesos*A)*A';

    l_adj=l+v;
    for i=1:n_angaz
        l_adj=l_adj*180/pi;
    end
        
    l_adj;

    Ql_adj=Ql-Qv;

    sigma_posteriori_2=v'*Pesos*v/gl;
        
    disp 'Teste da razao de variancias:'
    disp '                             '
 
    qui=icdf('chi2',1-nsign,gl);
        
    if (sigma_posteriori_2*gl/sigma_priori_2 < qui)
        disp 'Aceitar o ajustamento'
    else
        disp 'Rejeitar o ajustamento'
    end

% elipses de erro e/ou elipses de confianca    
    
end
    
% *************************************************************************
% ******************* OPCAO = 4 *******************************************
% *************************************************************************

if (opcao1 == 4)

    iteracao=1;

    clc
    disp '****************************************************************'
    disp '* Datum com minimizacao parcial do traco da matriz cofactor Qx *'
    disp '****************************************************************'    
    disp '                                                              '

    for j=1:n_mp;
       if (j<10)
           str=['        ',num2str(j),' - ',char(ponto(j))];
        else
            str=['       ',num2str(j),' - ',char(ponto(j))];        
        end
        disp (str)
        j=j+1;
    end

    disp '                                                          '

    num=input ('Indique o numero de pontos do datum: ');
    disp '(cada ponto contribui com 2 parametros: total parametros>d)'
    if (2*num < nulidade+1)
        disp '                                                      '
        disp 'Numero de pontos insuficiente...'
        %break;
    end

    disp '                                                          '

    for i=1:num
        datum1=input ('Indique um ponto seleccionado para datum: ');
        datum(i)=datum1;
    end
  
    for i=1:num
        corta(2*i-1)=2*datum(i)-1;
        corta(2*i)=2*datum(i);
    end
    
    if (n_diraz>0)
        R=zeros(2*n_mp+ne,2*num);  
        B=zeros(2*n_mp+ne,2*n_mp+ne);
        S=zeros(2*n_mp+ne,2*num);
    else
        R=zeros(2*n_mp,2*num);
        B=zeros(2*n_mp,2*n_mp);
        S=zeros(2*n_mp,2*num);
    end      
        
% inicio do processo iterativo

    for iteracao=1:maxit

        if (iteracao>1)
            % angaz
            for j=1:n_angaz;
                func=ang_az_lin(X(2*em_estacao(j)-1),X(2*em_estacao(j)),X(2*de_estacao(j)-1),X(2*de_estacao(j)),X(2*para_estacao(j)-1),X(2*para_estacao(j)));
                A(j,2*em_estacao(j)-1)=func(1,1);
                A(j,2*em_estacao(j))=func(2,1);
                A(j,2*de_estacao(j)-1)=func(3,1);
                A(j,2*de_estacao(j))=func(4,1);
                A(j,2*para_estacao(j)-1)=func(5,1);
                A(j,2*para_estacao(j))=func(6,1);
                j=j+1;
            end
            % disth
            for j=n_angaz+1:n_angaz+n_disth
                func=dist_hor_lin(X(2*de_distancia(j-n_angaz)-1),X(2*de_distancia(j-n_angaz)),X(2*para_distancia(j-n_angaz)-1),X(2*para_distancia(j-n_angaz)));
                A(j,2*de_distancia(j-n_angaz)-1)=func(1,1);
                A(j,2*de_distancia(j-n_angaz))=func(2,1);
                A(j,2*para_distancia(j-n_angaz)-1)=func(3,1);
                A(j,2*para_distancia(j-n_angaz))=func(4,1);
                j=j+1; 
            end
            % rumo
            for j=n_angaz+n_disth+1:n_angaz+n_disth+n_rumo
                func=rumo_lin(X(2*de_rumo(j-n_angaz-n_disth)-1),X(2*de_rumo(j-n_angaz-n_disth)),X(2*para_rumo(j-n_angaz-n_disth)-1),X(2*para_rumo(j-n_angaz-n_disth)));
                A(j,2*de_rumo(j-n_angaz-n_disth)-1)=func(1,1);
                A(j,2*de_rumo(j-n_angaz-n_disth))=func(2,1);
                A(j,2*para_rumo(j-n_angaz-n_disth)-1)=func(3,1);
                A(j,2*para_rumo(j-n_angaz-n_disth))=func(4,1);
                j=j+1;
            end
            % diraz   
            for j=n_angaz+n_disth+n_rumo+1:total_obs
                func=direccao_lin(X(2*em_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*em_direccao(j-n_angaz-n_disth-n_rumo)),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)-1),X(2*para_direccao(j-n_angaz-n_disth-n_rumo)));
                A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(1,1);
                A(j,2*em_direccao(j-n_angaz-n_disth-n_rumo))=func(2,1);
                A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo)-1)=func(3,1);
                A(j,2*para_direccao(j-n_angaz-n_disth-n_rumo))=func(4,1);
                j=j+1;
            end
            % colunas da matriz A que contem o coeficiente do R0 em cada estacao
            k=1;
            j1=1;
            for j=2:n_diraz
                if (em_direccao(j) == em_direccao(j-1))
                    A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1; 
                else
                    A(n_angaz+n_disth+n_rumo+j-1,2*n_mp+k)=1;
                    k=k+1;
                end 
                j=j+1;
            end
            if (n_diraz>0)
                A(n_angaz+n_disth+n_rumo+n_diraz,2*n_mp+k)=1;
            end
        end
       
% para eliminar uma coluna, utilizar o operador empty []

        A1=A;
                    
        A1(:,corta)=[];
  
        A1;
    
        A2=A(:,corta);
            
        A2;
     
% calculo da matriz L tal que A1*L=A2 (A1_1 e a inversa generalizada de A1)

        A1_1=pinv(A1);

        L=A1_1*A2;
                
% calculo da matriz S
        
        for i=1:2*num
            S(corta(i),i)=-1;
        end
        
        if (num == 3)
            for i=1:(2*datum(1)-2)
                for j=1:2*num
                    S(i,j)=L(i,j);
                end
            end
            for i=(2*datum(1)+1):(2*datum(2)-2)
                for j=1:2*num
                    S(i,j)=L(i-2,j);
                end
            end
            for i=(2*datum(2)+1):(2*datum(3)-2)
                for j=1:2*num
                    S(i,j)=L(i-4,j);
                end
            end
            for i=(2*datum(3)+1):2*n_mp
                for j=1:2*num
                    S(i,j)=L(i-6,j);
                end
            end
        end

        if (num == 4)
            for i=1:(2*datum(1)-2)
                for j=1:2*num
                    S(i,j)=L(i,j);
                end
            end
            for i=(2*datum(1)+1):(2*datum(2)-2)
                for j=1:2*num
                    S(i,j)=L(i-2,j);
                end
            end
            for i=(2*datum(2)+1):(2*datum(3)-2)
                for j=1:2*num
                    S(i,j)=L(i-4,j);
                end
            end
            for i=(2*datum(3)+1):(2*datum(4)-2)
                for j=1:2*num
                    S(i,j)=L(i-6,j);
                end
            end            
            for i=(2*datum(4)+1):2*n_mp
                for j=1:2*num
                    S(i,j)=L(i-8,j);
                end
            end
        end
        
        if (num == 5)
            for i=1:(2*datum(1)-2)
                for j=1:2*num
                    S(i,j)=L(i,j);
                end
            end
            for i=(2*datum(1)+1):(2*datum(2)-2)
                for j=1:2*num
                    S(i,j)=L(i-2,j);
                end
            end
            for i=(2*datum(2)+1):(2*datum(3)-2)
                for j=1:2*num
                    S(i,j)=L(i-4,j);
                end
            end
            for i=(2*datum(3)+1):(2*datum(4)-2)
                for j=1:2*num
                    S(i,j)=L(i-6,j);
                end
            end  
            for i=(2*datum(4)+1):(2*datum(5)-2)
                for j=1:2*num
                    S(i,j)=L(i-8,j);
                end
            end             
            for i=(2*datum(5)+1):2*n_mp
                for j=1:2*num
                    S(i,j)=L(i-10,j);
                end
            end
        end
 
        if (num == 6)
            for i=1:(2*datum(1)-2)
                for j=1:2*num
                    S(i,j)=L(i,j);
                end
            end
            for i=(2*datum(1)+1):(2*datum(2)-2)
                for j=1:2*num
                    S(i,j)=L(i-2,j);
                end
            end
            for i=(2*datum(2)+1):(2*datum(3)-2)
                for j=1:2*num
                    S(i,j)=L(i-4,j);
                end
            end
            for i=(2*datum(3)+1):(2*datum(4)-2)
                for j=1:2*num
                    S(i,j)=L(i-6,j);
                end
            end  
            for i=(2*datum(4)+1):(2*datum(5)-2)
                for j=1:2*num
                    S(i,j)=L(i-8,j);
                end
            end       
            for i=(2*datum(5)+1):(2*datum(6)-2)
                for j=1:2*num
                    S(i,j)=L(i-10,j);
                end
            end            
            for i=(2*datum(6)+1):2*n_mp
                for j=1:2*num
                    S(i,j)=L(i-12,j);
                end
            end
        end
        
        if (num == 7)
            for i=1:(2*datum(1)-2)
                for j=1:2*num
                    S(i,j)=L(i,j);
                end
            end
            for i=(2*datum(1)+1):(2*datum(2)-2)
                for j=1:2*num
                    S(i,j)=L(i-2,j);
                end
            end
            for i=(2*datum(2)+1):(2*datum(3)-2)
                for j=1:2*num
                    S(i,j)=L(i-4,j);
                end
            end
            for i=(2*datum(3)+1):(2*datum(4)-2)
                for j=1:2*num
                    S(i,j)=L(i-6,j);
                end
            end  
            for i=(2*datum(4)+1):(2*datum(5)-2)
                for j=1:2*num
                    S(i,j)=L(i-8,j);
                end
            end       
            for i=(2*datum(5)+1):(2*datum(6)-2)
                for j=1:2*num
                    S(i,j)=L(i-10,j);
                end
            end  
            for i=(2*datum(6)+1):(2*datum(7)-2)
                for j=1:2*num
                    S(i,j)=L(i-12,j);
                end
            end            
            for i=(2*datum(7)+1):2*n_mp
                for j=1:2*num
                    S(i,j)=L(i-14,j);
                end
            end
        end        
        
        if (num > 7)
            disp ' num demasiado grande'            
            break;
        end
        
        S;
  
% calculo da matriz B

        for i=1:2*num
            B(corta(i),corta(i))=-1;
        end     
        
        B;
        
% calculo da matriz R dos constrangimentos        
        
        R=B*S;

% calculo do vector fecho

        for j=1:n_disth
            dist_c(j)=dist_hor_cal(X(2*de_distancia(j)-1),X(2*de_distancia(j)),X(2*para_distancia(j)-1),X(2*para_distancia(j)));
            j=j+1;
        end
        
        if (n_disth>0)
            dist_c=dist_c';
        end
              
% transformacao de direccoes azimutais em rumos atraves do calculo do R0 em
% cada estacao; para o efeito, e necessario contar o numero de estacoes
% ocupadas e em cada estacao o numero de pontos visados; em  cada estacao
% calcula-se o rumo ponto estacao-ponto visado, tendo-se entao R0=REV-LEV

% calculo do R0 para cada leitura azimutal
    
        for j=1:n_diraz
            R0(j)=rumo(M(em_direccao(j)),P(em_direccao(j)),M(para_direccao(j)),P(para_direccao(j)))-diraz(j);
            if (R0(j)<0)
                R0(j)=R0(j)+2*pi;
            end
        end

% calculo do R0 medio e respectivo desvio padrao em cada estacao 

        sum_rumo=R0(1);
        sum_rumo2=sum_rumo^2;
        num_pontarias=1;
        j1=1;
        for j=2:n_diraz
            if (em_direccao(j) == em_direccao(j-1))
                sum_rumo=sum_rumo+R0(j);
                sum_rumo2=sum_rumo2+R0(j)^2;
                num_pontarias=num_pontarias+1;
            else
                j2=j-1;
                rumo_medio=sum_rumo/num_pontarias;        
                std_rumo_medio=sqrt((sum_rumo2-sum_rumo^2/(num_pontarias))/(num_pontarias-1));
                for k=j1:j2 
                    R0_medio(k)=rumo_medio;
                    sigma_R0_medio(k)=std_rumo_medio;
                    k=k+1;
                end
                sum_rumo=R0(j);
                sum_rumo2=sum_rumo^2;
                num_pontarias=1;
                j1=j;
            end 
            j=j+1;
        end

        j2=j-1;
        rumo_medio=sum_rumo/num_pontarias;        
        std_rumo_medio=sqrt((sum_rumo2-sum_rumo^2/(num_pontarias))/(num_pontarias-1));  
        for k=j1:j2 
            R0_medio(k)=rumo_medio;
            sigma_R0_medio(k)=std_rumo_medio;
            k=k+1;
        end

        R0=R0';
        R0_medio=R0_medio';
        sigma_R0_medio=sigma_R0_medio';

% calculo das direccoes azimutais "calculadas", dadas como DIRAZ=R-R0

        j1=1;
        for j=2:n_diraz
            if (em_direccao(j) == em_direccao(j-1))
                diraz_c(j-1)=rumo_cal(X(2*em_direccao(j-1)-1),X(2*em_direccao(j-1)),X(2*para_direccao(j-1)-1),X(2*para_direccao(j-1)))-R0_medio(j1);
                if (diraz_c(j-1)<0)
                    diraz_c(j-1)=diraz_c(j-1)+2*pi;
                end
            else
                diraz_c(j-1)=rumo_cal(X(2*em_direccao(j-1)-1),X(2*em_direccao(j-1)),X(2*para_direccao(j-1)-1),X(2*para_direccao(j-1)))-R0_medio(j1); 
                if (diraz_c(j-1)<0)
                    diraz_c(j-1)=diraz_c(j-1)+2*pi;
                end                   
                j1=j;
            end 
            j=j+1;
        end
        diraz_c(n_diraz)=rumo_cal(X(2*em_direccao(n_diraz)-1),X(2*em_direccao(n_diraz)),X(2*para_direccao(n_diraz)-1),X(2*para_direccao(n_diraz)))-R0_medio(j1); 
        if (diraz_c(n_diraz)<0)
            diraz_c(n_diraz)=diraz_c(n_diraz)+2*pi;
        end     
                   
        diraz_c=diraz_c';
        
% calculo dos rumos utilizando as coordenadas aproximadas

        for j=1:n_rumo
            rumo_c=rumo_cal(X(2*de_rumo(j)-1),X(2*de_rumo(j)),X(2*para_rumo(j)-1),X(2*para_rumo(j)));
            j=j+1;
        end         
                      
% calculo do vector W (vector fecho)

        for j=1:n_angaz
            W(j,1)=angaz_c(j)-angaz(j);
            if (W(j,1)<-pi);
                W(j,1)=W(j,1)+2*pi;
            end        
            if (W(j,1)>pi);
                W(j,1)=W(j,1)-2*pi;
            end
            j=j+1;
        end           
        for j=n_angaz+1:n_angaz+n_disth
            W(j,1)=dist_c(j-n_angaz)-disth(j-n_angaz);
            j=j+1;
        end
        for j=n_angaz+n_disth+1:n_angaz+n_disth+n_rumo
            W(j,1)=rumo(j-n_angaz-n_disth)-rumo_c(j-n_angaz-n_disth);
            j=j+1;
        end        
        for j=1:n_diraz
            W(n_disth+j,1)=diraz_c(j)-diraz(j);
            if (W(n_disth+j,1)<-pi);
                W(n_disth+j,1)=W(n_disth+j,1)+2*pi;
            end        
            if (W(n_disth+j,1)>pi);
                W(n_disth+j,1)=W(n_disth+j,1)-2*pi;
            end
            j=j+1;
        end
    
        W;     
        
% calculo do vector dos parametros

        B1=inv(N+R*R')*N*inv(N+R*R');

        dx=-B1*A'*Pesos*W;
           
% calculo das novas coordenadas aproximadas

        if (max(abs(dx)) > conv)
            for i=1:2*n_mp; 
                X(i,1)=X(i,1)+dx(i,1);
                i=i+1;
            end
        else
            disp ('O processo iterativo convergiu.')
            disp ('                                                    ')
            break;
        end
        
        if (iteracao > maxit)
            disp ('O processo iterativo nao convergiu ...')
            stop       
        end
            
        iteracao=iteracao+1;
    
    end
    
    disp 'Vector das coordenadas ajustadas:'
    X
    
    QX=inv(N+R*R')*N*inv(N+R*R');

    tracoQx=0;
    
    for i=1:2*n_mp
        tracoQx=tracoQx+QX(i,i);
    end

    disp 'Traco da matriz cofactor das coordenadas:'
    tracoQx

    disp 'Vector dos residuos das observacoes:'
    v=A*dx+W

    Qv=Ql-A*inv(A'*Pesos*A)*A';

    l_adj=l+v;
    for i=1:n_angaz
        l_adj=l_adj*180/pi;
    end
        
    l_adj;

    Ql_adj=Ql-Qv;

    sigma_posteriori_2=v'*Pesos*v/gl;
        
    disp 'Teste da razao de variancias:'
    disp '                             '
 
    qui=icdf('chi2',1-nsign,gl);
        
    if (sigma_posteriori_2*gl/sigma_priori_2 < qui)
        disp 'Aceitar o ajustamento'
    else
        disp 'Rejeitar o ajustamento'
    end

% elipses de erro e/ou elipses de confianca        
    
end

% *************************************************************************
% ******************* OPCAO = 5 *******************************************
% *************************************************************************

if (opcao1 == 5)
end

% elipses de erro

fid=fopen('elipses.txt','w');

for i=1:n_mp
    a=sqrt(((QX(2*i-1,2*i-1)+QX(2*i,2*i)+sqrt((QX(2*i-1,2*i-1)+QX(2*i,2*i))^2-4*(QX(2*i-1,2*i-1)*QX(2*i,2*i)-QX(2*i-1,2*i)^2))))/2);
    b=sqrt(((QX(2*i-1,2*i-1)+QX(2*i,2*i)-sqrt((QX(2*i-1,2*i-1)+QX(2*i,2*i))^2-4*(QX(2*i-1,2*i-1)*QX(2*i,2*i)-QX(2*i-1,2*i)^2))))/2);
% teta = angulo entre o eixo Y e o semi-eixo maior da elipse, no sentido horario        
    if (QX(2*i-1,2*i-1) == QX(2*i,2*i))
        teta=0; %circunferencia
    elseif (QX(2*i-1,2*i) == 0)
        if (QX(2*i-1,2*i) > QX(2*i,2*i))
            teta=pi/2;
        else
            teta=0;
        end
    else
        teta=0.5*atan2(2*QX(2*i-1,2*i),QX(2*i,2*i)-QX(2*i-1,2*i-1));
    end
    if (teta < 0)
        teta=teta+pi;
    end
    if (a ~= 0 & b ~= 0)
        %fprintf(fid,'%11.4f %11.4f %11.4f %11.4f %11.4f\n',X(2*i-1),X(2*i),a*esc*conf,b*esc*conf,teta*180/pi);
    end
end    
    
fclose(fid);
