К СКН

[eddy_em]

В прошлой теме по поводу СКН БТА натолкнули меня на простейшую мысль решать систему уравнений совместно.
Попробовал: вполне себе решается, и решение даже не очень сильно восприимчиво к шумам. Вот, например, беру реальные данные и добавляю разный шум:

	Noise = 0.000000, sigmaA=0.000000, sigmaZ=0.000000
PCS coefficients: K0..K7: -92.747067, 12.012553, 32.525176, 2.286829, 3.233337, -97.918763, 13.777403, 12.757328
Differencies: min=-12.6, max=27.9, std=5.8
	Noise = 0.500000, sigmaA=0.281029, sigmaZ=0.271943
PCS coefficients: K0..K7: -93.158939, 11.431714, 33.213620, 2.301323, 3.200591, -98.126947, 13.850065, 12.886362
Differencies: min=-12.8, max=27.7, std=5.8
	Noise = 1.000000, sigmaA=0.522671, sigmaZ=0.580003
PCS coefficients: K0..K7: -93.227872, 11.181730, 33.431938, 2.253908, 3.217875, -97.966544, 14.024323, 12.616600
Differencies: min=-11.8, max=27.7, std=5.7
	Noise = 3.000000, sigmaA=1.539450, sigmaZ=1.731406
PCS coefficients: K0..K7: -89.141946, 16.210061, 27.572613, 2.206247, 3.432656, -107.467989, 21.205184, 20.001838
Differencies: min=-14.1, max=29.7, std=6.1
	Noise = 5.000000, sigmaA=2.955128, sigmaZ=2.880220
PCS coefficients: K0..K7: -90.050438, 17.520024, 26.707310, 2.372867, 2.963566, -87.841441, 7.462785, 4.779566
Differencies: min=-13.3, max=30.8, std=6.7

Здесь Noise — размах равномерно распределенного аддитивного шума, sigmaX — RMS отклонения зашумленных данных от реальных, Differencies — остаточные невязки.
Данные собраны были по всего-то тридцати точкам, поэтому результаты не очень.

function PCS = getPCScoeff(tabname, imprefix="", nheader=8)
#{
PCS = getPCScoeff(tabname, imprefix, nheader)

Calculate PCS coefficients & plot graphs

Parameters:
    tabname - filename with table of dA/dZ
    imprefix - prefix for generated graphics ("" by default)
    nheader - N lines of header file in `tabname` (8 by default)

Pointing Correction System coefficients:
 dA = K0 + K1/tg(Z) + K2/sin(Z) - K3*sin(A)/tg(Z) + K4 *cos(delta)*cos(P)/sin(Z)
 dZ = K5 + K6*sin(Z) + K7*cos(Z) + K3*cos(A) + K4*cos(phi)*sin(A)

 K0 = A0 - azimuth zero; K1 = L - horiz axe inclination; K2 = k - collimation error;
 K3 = F - lattitude error of vert. axe; K4 = dS - time error
 K5 = Z0 - zenith zero; K6 = d - tube bend; K7 = d1 - cos. tube bend

 phi = 43.6535278 - lattitude
 t = LST - Alpha  - hour angle
 P=atan(sin(t)/(tan(phi)*cos(Del)-sin(Del)*cos(t)))  - parallax angle
#}
    f = fopen(tabname);
    ALL = textscan(f, "|%f:%f:%f %f:%f:%f|%f %f|%f %f|%f %f|%f:%f:%f|", "HeaderLines", 8);
    fclose(f);
    %     1    2    3  4    5    6    7     8      9 10 11 12  13  14  15
    %    [Ald Alm Als Deld Delm Dels dAl_S dDel_S dA dZ  A  Z  STh STm STs ]
    A = ALL{11}*pi/180; % all angles here will be in radians
    Z = ALL{12}*pi/180;
    Al = pi*(ALL{1}+ALL{2}/60+ALL{3}/3600)/12; % right accession from hours to rad
    Delsig = ALL{4}./abs(ALL{4});   % declination sign
    Del = pi*Delsig.*(abs(ALL{4})+ALL{5}/60+ALL{6}/3600)/180; % declination from degrees to rad
    phi = 43.6535278 * pi / 180; % lattitude
    t = pi*(ALL{13}+ALL{14}/60+ALL{15}/3600)/12 - Al; % hour angle
    P = atan(sin(t)./(tan(phi).*cos(Del)-sin(Del).*cos(t))); % parallax angle
    cosZ = cos(Z);
    sinZ = sin(Z);
    cosA = cos(A);
    sinA = sin(A);
    tgZ  = tan(Z);
    dA = ALL{9}; % dA/dZ in arcsec
    dZ = ALL{10};
    halfsz = size(dA); 
    ZRS = zeros(halfsz);
    % MATRIX * [K0 K1 K2 K3 K4 K5 K6 K7 ]' = [dA dZ]'
    Matrix = [ ...
        [ones(halfsz); ZRS]         ... % K0 - dA
        [1./tgZ; ZRS]               ... % K1/tg(Z) - dA
        [1./sinZ; ZRS]              ... % K2/sin(Z) - dA
        [sinA./tgZ; cosA]           ... % K3*sin(A)/tg(Z) - dA, K3*cos(A) - dZ
        [cos(Del).*cos(P)./sinZ; cos(phi)*sinA] ... % K4 *cos(delta)*cos(P)/sin(Z) - dA, K4*cos(phi)*sin(A) - dZ
        [ZRS; ones(halfsz)]         ... % K5 - dZ
        [ZRS; sinZ]                 ... % K6*sin(Z) - dZ
        [ZRS; cosZ]                 ... % K7*cos(Z) - dZ
    ];
    PCS = Matrix \ [dA; dZ];
    for N = [0 0.5 1 3 5]; 
        ddA = dA + 2*N*(rand(halfsz)-0.5);
        ddZ = dZ + 2*N*(rand(halfsz)-0.5);
        printf("\tNoise = %f, sigmaA=%f, sigmaZ=%f\n", N, std(ddA-dA), std(ddZ-dZ));
        getK(Matrix, ddA, ddZ);
    endfor
    if(isempty(imprefix)) return; endif
    diffs = Matrix * PCS - [dA; dZ];
    L = length(diffs);
    ddA = diffs(1:L/2); ddZ = diffs(L/2+1:L);
    fg = figure;
    plot(A*180/pi, [ddA ddZ], 'o');
    legend("ddA", "ddZ");
    xlabel("A, degr"); ylabel("Remaining error: real-model");
    plotgr(sprintf("%s_%s", imprefix, "diff_vs_A"), fg);
    fg = figure;
    plot(Z*180/pi, [ddA ddZ], 'o');
    legend("ddA", "ddZ");
    xlabel("Z, degr"); ylabel("Remaining error: real-model");
    plotgr(sprintf("%s_%s", imprefix, "diff_vs_Z"), fg);
    #}
endfunction

function K = getK(M, dA, dZ)
    K = M \ [dA; dZ];
    printf("PCS coefficients: K0..K7: %f, %f, %f, %f, %f, %f, %f, %f\n", ...
        K(1), K(2), K(3), K(4), K(5), K(6), K(7), K(8));
    diffs = M * K - [dA; dZ];
    printf("Differencies: min=%.1f, max=%.1f, std=%.1f\n", min(diffs), max(diffs), std(diffs));
endfunction

function plotgr(nm, fg)
    print(fg, '-dpdf', sprintf("%s.pdf", nm));
    print(fg, '-dpng', sprintf("%s.png", nm));
endfunction

Он работает со стандартным файлом отчета, поэтому такая дичь в предварительной подготовке. Зато визуально более-менее понятно.
Никак не могу найти более полных наборов данных — чтобы хотя бы точек 80 было. Нашел данные за этот год — 86 точек, надо будет обработать. Заодно соберу в одном месте все нужные скрипты, а то каждый раз приходится искать, где же версия посвежей! Да и на гит это надо закинуть.