#!/bin/bash

###############################################################################

ct="(cos(tau*pi()/180.0))"
st="(sin(tau*pi()/180.0))"
cd="(cos(dec*pi()/180.0))"
sd="(sin(dec*pi()/180.0))"

ct0="(cos(tt*pi()/180.0))"
st0="(sin(tt*pi()/180.0))"
cd0="(cos(dd*pi()/180.0))"
sd0="(sin(dd*pi()/180.0))"

fx="           +b0*$sd    -c0*$cd*$st+d0*$sd*$st+e0*$sd*$ct-f0*$st"
fy="-a0*$sd               +c0*$cd*$ct-d0*$sd*$ct+e0*$sd*$st+f0*$ct"
fz="+a0*$cd*$st-b0*$cd*$ct                      -e0*$cd"

# $cd*$st, $cd*$ct, $sd

y1="($cd*$st)"
y2="($sd)"
y3="($cd*$ct)"

fx="$fx+b1*$y1*$sd-c1*$y1*$cd*$st+d1*$y1*$sd*$st+e1*$y1*$sd*$ct-f1*$y1*$st"
fy="$fy-a1*$y1*$sd+c1*$y1*$cd*$ct-d1*$y1*$sd*$ct+e1*$y1*$sd*$st+f1*$y1*$ct"
fz="$fz+a1*$y1*$cd*$st-b1*$y1*$cd*$ct-e1*$y1*$cd"

fx="$fx+b2*$y2*$sd-c2*$y2*$cd*$st+d2*$y2*$sd*$st+e2*$y2*$sd*$ct-0*$y2*$st"
fy="$fy-a2*$y2*$sd+c2*$y2*$cd*$ct-d2*$y2*$sd*$ct+e2*$y2*$sd*$st+0*$y2*$ct"
fz="$fz+a2*$y2*$cd*$st-b2*$y2*$cd*$ct-e2*$y2*$cd"

fx="$fx+b3*$y3*$sd-c3*$y3*$cd*$st+d3*$y3*$sd*$st+e3*$y3*$sd*$ct-f3*$y3*$st"
fy="$fy- 0*$y3*$sd+c3*$y3*$cd*$ct-d3*$y3*$sd*$ct+e3*$y3*$sd*$st+f3*$y3*$ct"
fz="$fz+ 0*$y3*$cd*$st-b3*$y3*$cd*$ct-e3*$y3*$cd"

yx="($cd0*$ct0-$cd*$ct)"
yy="($cd0*$st0-$cd*$st)"
yz="($sd0-$sd)"

cdef=tau:1,dec:2,tt:3,dd:4


# change 0 to 1 in ( ... || 0 ) below in order to include the points
# after pole crossing (i.e. having an encoder declination larger than 90).

# first, check the residual of the pointing without any pointing model:
if false; then
	cat	pointing.txt | \
	awk '{ if ( $2 <= 90 || 0 ) print; }' | \
	lfit	-c $cdef -f "$yx^2+$yy^2+$yz^2" | \
	awk '{ sum+=$1;n++; } END { printf("%10.6f\n",sqrt(sum/n)); }'
fi

cat	pointing.txt | \
awk '{ if ( $2 <= 90 || 1 ) print; }' | \
lfit	-v a0,b0,c0,d0,e0,f0,a1,b1,c1,d1,e1,f1,a2,b2,c2,d2,e2,b3,c3,d3,e3,f3 \
	-F a0=%9.6f,b0=%9.6f,c0=%9.6f,d0=%9.6f,e0=%9.6f,f0=%9.6f \
	-F a1=%9.6f,b1=%9.6f,c1=%9.6f,d1=%9.6f,e1=%9.6f,f1=%9.6f \
	-F a2=%9.6f,b2=%9.6f,c2=%9.6f,d2=%9.6f,e2=%9.6f \
	-F b3=%9.6f,c3=%9.6f,d3=%9.6f,e3=%9.6f,f3=%9.6f \
	-ix - -cx $cdef -fx "$fx" -yx "$yx" -ex 0.000013 \
	-iy - -cy $cdef -fy "$fy" -yy "$yy" -ey 0.000013 \
	-iz - -cz $cdef -fz "$fz" -yz "$yz" -ez 0.000013 \
	--err --residual | tee pointing.fit

if true; then
	A=($(cat pointing.fit))
	cat	pointing.txt | \
	lfit	-v a0=${A[0]},b0=${A[1]},c0=${A[2]},d0=${A[3]},e0=${A[4]},f0=${A[5]},a1=${A[6]},b1=${A[7]},c1=${A[8]},d1=${A[9]},e1=${A[10]},f1=${A[11]},a2=${A[12]},b2=${A[13]},c2=${A[14]},d2=${A[15]},e2=${A[16]},b3=${A[17]},c3=${A[18]},d3=${A[19]},e3=${A[20]},f3=${A[21]} \
		-c $cdef \
		-f "tau,dec,sqrt(($fx-$yx)^2+($fy-$yy)^2+($fz-$yz)^2),$fx-$yx,$fy-$yy,$fz-$yz" \
		-F %10.5f,%10.5f,%10.6f,%10.6f,%10.6f,%10.6f \
		-o pointing.res
fi

###############################################################################

