Astronomy: interpolation

[salvomic]

hi Marcel and every friends,
about Astronomy it's very important to have some tools to get interpolation (in part used in our Astro Lab and Effemeridi for Prime).
I think this is the base for an ancillary program for almost every task.
Please, help to improve and doing the necessary controls to avoid logic errors and crashes too...

The tools are:
interpolation of 3 or 5 values
finding the "n" interpolation factor
finding the extremum (minimum, maximum of a list to interpolate)
finding the zero (3 or 5 values)
"halves" (4 values)
get a value with Lagrange interpolation.

Typical example: interpolation({1.9556788, 1.97345409, 1.965049586}, 0.18), nzero3({1.9556788, 1.97345409, 1.965049586}) ...
Very useful for RA (right ascension) and Declination (or lambda, beta: longitude, latitude) for an aster or planet or to find conjunction of planets, least angular distance and so on...
Formulae by Jean Meeus, Astronomical Algorithms, very precise, even more with the power of the Prime!

Salvo

Code:


// Ancillary program for Astro Lab 4 (Marcel Pelletier) and Effemeridi - Astrolabio (Salvo Micciché)
// lista = a list like {data1, data2, data3...}
// n = interpolation factor

EXPORT interpolation3(lista,n)
BEGIN
  LOCAL a,b,c,R,dif1,dif2;
  dif1:=ΔLIST(lista);
  dif2:=ΔLIST(dif1);
  a:=dif1(1);
  b:=dif1(2);
  c:=dif2(1);
  R:=lista(2)+(n/2)*(a+b+n*c);
  RETURN R;
END;

EXPORT interpolation5(lista,n)
BEGIN
  LOCAL a,b,c,R,dif1,dif2, dif3, dif4;
  LOCAL d,e,f,g,h,j,k;
  dif1:=ΔLIST(lista);
  dif2:=ΔLIST(dif1);
  dif3:=ΔLIST(dif2);
  dif4:=ΔLIST(dif3);
  a:=dif1(1);
  b:=dif1(2);
  c:=dif1(3);
  d:=dif1(4);
  e:=dif2(1);
  f:=dif2(2);
  g:=dif2(3);
  h:=dif3(1);
  j:=dif3(2);
  k:=dif4(1);
  R:=lista(3)+(n/2)*(b+c)+(n^2/2)*f+n*((n^2-1)/12)*(h+j)+(n^2*(n^2-1)/24)*k;
  RETURN R;
END;

EXPORT extremum3(lista)
BEGIN
  LOCAL a,b,c,R,dif1,dif2, ym, nm;
  dif1:=ΔLIST(lista);
  dif2:=ΔLIST(dif1);
  a:=dif1(1);
  b:=dif1(2);
  c:=dif2(1);
  ym:=lista(2)-((a+b)^2)/(8*c);
  nm:= -(a+b)/(2*c);
  RETURN ({ym, nm});
END;

EXPORT extremum5(lista)
BEGIN
  LOCAL a,b,c,R,dif1,dif2, dif3, dif4;
  LOCAL d,e,f,g,h,j,k, nm;
  LOCAL jj, temp;
  dif1:=ΔLIST(lista);
  dif2:=ΔLIST(dif1);
  dif3:=ΔLIST(dif2);
  dif4:=ΔLIST(dif3);
  a:=dif1(1);
  b:=dif1(2);
  c:=dif1(3);
  d:=dif1(4);
  e:=dif2(1);
  f:=dif2(2);
  g:=dif2(3);
  h:=dif3(1);
  j:=dif3(2);
  k:=dif4(1);
  jj:= 0; nm:=0; temp:= 0;
  WHILE 1 DO
  temp:= nm;
  nm:= (6*b+6*c-h-j+3*jj^2*(h+j)+2*jj^3*k)/(k-12*f);
    IF nm==temp THEN BREAK; END;
  jj:= nm;
  END;
  RETURN nm;
END;

EXPORT nzero3(lista)
BEGIN
  LOCAL a,b,c,R,dif1,dif2, n0, j;
  LOCAL temp;
  dif1:=ΔLIST(lista);
  dif2:=ΔLIST(dif1);
  a:=dif1(1);
  b:=dif1(2);
  c:=dif2(1);
  j:= 0; n0:=0; temp:= 0;
  WHILE 1 DO
  temp:= n0;
  n0:= -2*lista(2)/(a+b+c*j);
    IF n0==temp THEN BREAK; END;
  j:= n0;
  END;
  RETURN n0;
END;

EXPORT nzero5(lista)
BEGIN
  LOCAL a,b,c,R,dif1,dif2, dif3, dif4;
  LOCAL d,e,f,g,h,j,k, n0;
  LOCAL jj, temp;
  dif1:=ΔLIST(lista);
  dif2:=ΔLIST(dif1);
  dif3:=ΔLIST(dif2);
  dif4:=ΔLIST(dif3);
  a:=dif1(1);
  b:=dif1(2);
  c:=dif1(3);
  d:=dif1(4);
  e:=dif2(1);
  f:=dif2(2);
  g:=dif2(3);
  h:=dif3(1);
  j:=dif3(2);
  k:=dif4(1);
  jj:= 0; n0:=0; temp:= 0;
  WHILE 1 DO
  temp:= n0;
  n0:= (-24*lista(3)+jj^2*(k-12*f)-2*jj^3*(h+j)-jj^4*k)/(2*(6*b+6*c-h-j));
    IF n0==temp THEN BREAK; END;
  jj:= n0;
  END;
  RETURN n0;
END;

EXPORT halves4(lista)
BEGIN
  LOCAL y;
  y:= (9*(lista(2)+lista(3))-lista(1)-lista(4))/16;
  RETURN y;
END;

EXPORT lagrangeValue(lista_x, lista_y, xValue)
BEGIN
  LOCAL y, temp;
  temp:= lagrange(lista_x, lista_y);
  y:= EVAL(EXPR(temp+"|X="+xValue));
  RETURN y;
END;