From a860d08c11a5ed1ff3174d3aa44a49c4920e00db Mon Sep 17 00:00:00 2001
From: Sarod Yatawatta <sarod@dop373.nfra.nl>
Date: Thu, 18 Apr 2024 00:31:22 +0200
Subject: [PATCH 01/15] add cspice

---
 CMakeLists.txt        | 14 ++++++++++++++
 src/MS/CMakeLists.txt |  1 +
 2 files changed, 15 insertions(+)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index ea7269f..73e3e47 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -123,6 +123,20 @@ pkg_check_modules(GLIB_PKG glib-2.0)
 if (GLIB_PKG_FOUND)
     include_directories(${GLIB_PKG_INCLUDE_DIRS})
 endif()
+
+#CSPICE find if -DCSPICE_PREFIX=... given
+set(CSPICE_PREFIX
+  ""
+  CACHE FILEPATH "Path to CSPICE install root")
+if(NOT "${CSPICE_PREFIX}" STREQUAL "")
+  set(ENV{PKG_CONFIG_PATH} "${CSPICE_PREFIX}/lib/pkgconfig")
+  pkg_search_module(CSPICE libcspice)
+endif()
+if(CSPICE_FOUND)
+  message(STATUS "Found CSPICE: ${CSPICE_INCLUDE_DIRS}")
+  include_directories(${CSPICE_INCLUDE_DIRS})
+  add_definitions(-DHAVE_CSPICE)
+endif()
 #--------------------------------------- check environment variables
 
 # Check pre-defined environment variables
diff --git a/src/MS/CMakeLists.txt b/src/MS/CMakeLists.txt
index 789e7df..22b724e 100644
--- a/src/MS/CMakeLists.txt
+++ b/src/MS/CMakeLists.txt
@@ -40,6 +40,7 @@ if(HAVE_CUDA)
         ${LAPACK_LIBRARIES}
         ${WCSLIB_LIBRARIES}
         ${GLIB_PKG_LIBRARIES}
+        $<$<1:${CSPICE_FOUND}>:${CSPICE_LINK_LIBRARIES}>
         -lpthread
         -lm
         ${CUDA_CUBLAS_LIBRARIES}

From 8644b6c266c64d7e0a641bd4d2f914448b2e5686 Mon Sep 17 00:00:00 2001
From: Sarod Yatawatta <sarod@dop373.nfra.nl>
Date: Thu, 18 Apr 2024 13:19:51 +0200
Subject: [PATCH 02/15] add uvwriter

---
 CMakeLists.txt              |  9 -----
 src/CMakeLists.txt          |  7 ++++
 src/MS/CMakeLists.txt       |  3 +-
 src/uvwriter/CMakeLists.txt | 68 +++++++++++++++++++++++++++++++++++++
 src/uvwriter/uvwriter.cpp   | 44 ++++++++++++++++++++++++
 5 files changed, 121 insertions(+), 10 deletions(-)
 create mode 100644 src/uvwriter/CMakeLists.txt
 create mode 100644 src/uvwriter/uvwriter.cpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 73e3e47..c8df1d6 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -195,8 +195,6 @@ message(STATUS "\n############################\n#   Configuration summary\n#####
 message (STATUS "CMAKE_SYSTEM .............. = ${CMAKE_SYSTEM}")
 message (STATUS "CMAKE_INSTALL_PREFIX ...... = ${CMAKE_INSTALL_PREFIX}")
 message (STATUS "CMAKE_BUILD_TYPE .......... = ${CMAKE_BUILD_TYPE}")
-# message (STATUS "BUILD_SHARED_LIBS ......... = ${BUILD_SHARED_LIBS}")
-# message (STATUS "USE_STACKTRACE ............ = ${USE_STACKTRACE}")
 message (STATUS "CMAKE_CXX_COMPILER ........ = ${CMAKE_CXX_COMPILER}")
 message (STATUS "CMAKE_CXX_FLAGS_RELEASE ... = ${CMAKE_CXX_FLAGS_RELEASE}")
 
@@ -238,11 +236,8 @@ message (STATUS "OpenBLAS_LIB .............. = ${OpenBLAS_LIB}")
 message (STATUS "GLIB_PKG_INCLUDE_DIRS...... = ${GLIB_PKG_INCLUDE_DIRS}")
 message (STATUS "GLIB_PKG_LIBRARIES......... = ${GLIB_PKG_LIBRARIES}")
 
-#message (STATUS "LAPACK-INC ................ = ${LAPACK_INCLUDE_DIR}")
 message (STATUS "LAPACK_LIBRARIES........... = ${LAPACK_LIBRARIES}")
 
-# message (STATUS "GFORTRAN-INC  ............. = ${GFORTRAN_INCLUDE_DIR}")
-# message (STATUS "GFORTRAN-LIBS ............. = ${LIBGFORTRAN_LIBRARIES}")
 message (STATUS "CMAKE_Fortran_COMPILER..... = ${CMAKE_Fortran_COMPILER}")
 
 message (STATUS "CFITSIO_ROOT_DIR........... = ${CFITSIO_ROOT_DIR}")
@@ -252,10 +247,6 @@ message (STATUS "CFITSIO_LIB................ = ${CFITSIO_LIB}")
 message (STATUS "WCSLIB_INCLUDE_DIRS........ = ${WCSLIB_INCLUDE_DIRS}")
 message (STATUS "WCSLIB_LIBRARIES .......... = ${WCSLIB_LIBRARIES}")
 
-# message (STATUS "HDF5_INCLUDE_DIR........... = ${HDF5_INCLUDE_DIRS}")
-# message (STATUS "HDF5_LIBRARIES............. = ${HDF5_LIBRARIES}")
-
-
 #--------------------------------------- include directories
 add_subdirectory(src)
 
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index e27f7fc..0dbdc9e 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -28,3 +28,10 @@ if(NOT LIB_ONLY)
    message(STATUS "\n WCSLIB not found, so skipping buildsky/restore")
   endif()
 endif()
+
+
+if(NOT LIB_ONLY)
+  if(CSPICE_FOUND)
+   add_subdirectory(uvwriter)
+  endif()
+endif()
diff --git a/src/MS/CMakeLists.txt b/src/MS/CMakeLists.txt
index 22b724e..039f6e1 100644
--- a/src/MS/CMakeLists.txt
+++ b/src/MS/CMakeLists.txt
@@ -40,7 +40,7 @@ if(HAVE_CUDA)
         ${LAPACK_LIBRARIES}
         ${WCSLIB_LIBRARIES}
         ${GLIB_PKG_LIBRARIES}
-        $<$<1:${CSPICE_FOUND}>:${CSPICE_LINK_LIBRARIES}>
+        $<$<BOOL:${CSPICE_FOUND}>:${CSPICE_LINK_LIBRARIES}>
         -lpthread
         -lm
         ${CUDA_CUBLAS_LIBRARIES}
@@ -66,6 +66,7 @@ else()
         ${LAPACK_LIBRARIES}
         ${WCSLIB_LIBRARIES}
         ${GLIB_PKG_LIBRARIES}
+        $<$<BOOL:${CSPICE_FOUND}>:${CSPICE_LINK_LIBRARIES}>
         -lpthread
         -lm
     )
diff --git a/src/uvwriter/CMakeLists.txt b/src/uvwriter/CMakeLists.txt
new file mode 100644
index 0000000..5949c03
--- /dev/null
+++ b/src/uvwriter/CMakeLists.txt
@@ -0,0 +1,68 @@
+include_directories(${CASACORE_INCLUDE_DIR})
+include_directories(${CASACORE_INCLUDE_DIR}/casacore)
+
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/../lib/Dirac)
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/../lib/Radio)
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/../lib/Radio/reserve)
+include_directories(${CMAKE_CURRENT_SOURCE_DIR})
+
+if(HAVE_CUDA)
+    CUDA_INCLUDE_DIRECTORIES(${CUDA_INCLUDE_DIRS})
+    CUDA_INCLUDE_DIRECTORIES(${CASACORE_INCLUDE_DIR})
+    CUDA_INCLUDE_DIRECTORIES(${CASACORE_INCLUDE_DIR}/casacore)
+    CUDA_INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR}/../lib/Dirac)
+    CUDA_INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR}/../lib/Radio)
+    CUDA_INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR}/../lib/Radio/reserve)
+    CUDA_INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR})
+endif()
+
+link_directories(${LIBRARY_OUTPUT_PATH})
+link_directories(${CMAKE_CURRENT_SOURCE_DIR}/../lib/Dirac)
+link_directories(${CMAKE_CURRENT_SOURCE_DIR}/../lib/Radio)
+
+FILE(GLOB SRCFILES *.cpp)
+
+# add build dependencies
+target_link_libraries(${target} dirac-radio)
+target_link_libraries(${target} dirac)
+
+
+if (HAVE_CUDA)
+  cuda_add_executable(uvwriter ${SRCFILES})
+  add_dependencies(uvwriter dirac-radio dirac)
+  target_link_libraries(uvwriter
+        -ldirac-radio
+        -ldirac
+        ${CASACORE_LIBRARIES}
+        ${CFITSIO_LIB}
+        ${OpenBLAS_LIB}
+        ${LAPACK_LIBRARIES}
+        ${WCSLIB_LIBRARIES}
+        ${GLIB_PKG_LIBRARIES}
+        $<$<1:${CSPICE_FOUND}>:${CSPICE_LINK_LIBRARIES}>
+        -lpthread
+        -lm
+        ${CUDA_CUBLAS_LIBRARIES}
+        ${CUDA_CUFFT_LIBRARIES}
+        ${CUDA_cusolver_LIBRARIES}
+        ${CUDA_cudadevrt_LIBRARIES}
+        ${NVML_LIB_PATH}
+  )
+else()
+  add_executable(uvwriter ${SRCFILES})
+  add_dependencies(uvwriter dirac-radio dirac)
+  target_link_libraries(uvwriter
+        -ldirac-radio
+        -ldirac
+        ${CASACORE_LIBRARIES}
+        ${CFITSIO_LIB}
+        ${OpenBLAS_LIB}
+        ${LAPACK_LIBRARIES}
+        ${WCSLIB_LIBRARIES}
+        ${GLIB_PKG_LIBRARIES}
+        $<$<1:${CSPICE_FOUND}>:${CSPICE_LINK_LIBRARIES}>
+        -lpthread
+        -lm
+  )
+endif()
+install(TARGETS uvwriter DESTINATION bin)
diff --git a/src/uvwriter/uvwriter.cpp b/src/uvwriter/uvwriter.cpp
new file mode 100644
index 0000000..a33616a
--- /dev/null
+++ b/src/uvwriter/uvwriter.cpp
@@ -0,0 +1,44 @@
+/*
+ *
+ Copyright (C) 2024 Sarod Yatawatta <sarod@users.sf.net>
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ $Id$
+ */
+
+#include <unistd.h>
+#include <stdio.h>
+#include <iostream>
+#include <casacore/ms/MeasurementSets/MSIter.h>
+#include <casacore/tables/Tables/Table.h>
+#include <casacore/tables/Tables/TableVector.h>
+#include <casacore/tables/Tables/TableRecord.h>
+#include <casacore/tables/Tables/TableColumn.h>
+#include <casacore/tables/Tables/TableIter.h>
+#include <casacore/tables/Tables/ScalarColumn.h>
+#include <casacore/tables/Tables/ArrayColumn.h>
+#include <casacore/casa/Arrays/Array.h>
+#include <casacore/casa/Arrays/Cube.h>
+#include <fstream>
+#include <math.h>
+#include <complex>
+#include <SpiceUsr.h>
+#include "Dirac_radio.h"
+
+
+int 
+main(int argc, char **argv) {
+
+  return 0;
+}

From a862b43bdcaa82bd5468ceec8bf0cecaa9cfa1cc Mon Sep 17 00:00:00 2001
From: Sarod Yatawatta <sarod@dop373.nfra.nl>
Date: Sat, 20 Apr 2024 23:07:50 +0200
Subject: [PATCH 03/15] add writer

---
 src/uvwriter/uvwriter.cpp | 210 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 210 insertions(+)

diff --git a/src/uvwriter/uvwriter.cpp b/src/uvwriter/uvwriter.cpp
index a33616a..b45cbe5 100644
--- a/src/uvwriter/uvwriter.cpp
+++ b/src/uvwriter/uvwriter.cpp
@@ -19,6 +19,7 @@
 
 #include <unistd.h>
 #include <stdio.h>
+#include <stdlib.h>
 #include <iostream>
 #include <casacore/ms/MeasurementSets/MSIter.h>
 #include <casacore/tables/Tables/Table.h>
@@ -36,9 +37,218 @@
 #include <SpiceUsr.h>
 #include "Dirac_radio.h"
 
+//#define DEBUG
+
+void
+cspice_load_kernels(void) {
+  /* get env CSPICE_KERNEL_PATH */
+  char* cspice_path = getenv("CSPICE_KERNEL_PATH");
+  if (cspice_path) {
+    const char *kname="/pck00011.tpc\0";
+    char *fullname=(char*)calloc((size_t)strlen((char*)cspice_path)+1+strlen((char*)kname),sizeof(char));
+    if (fullname == 0 ) {
+         fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
+         exit(1);
+    }
+    strcpy(fullname,cspice_path);
+    strcpy((char*)&(fullname[strlen(cspice_path)]),kname);
+    printf("loading %s\n",fullname);
+    furnsh_c(fullname);
+    free(fullname);
+
+    kname="/naif0012.tls\0";
+    fullname=(char*)calloc((size_t)strlen((char*)cspice_path)+1+strlen((char*)kname),sizeof(char));
+    if (fullname == 0 ) {
+         fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
+         exit(1);
+    }
+    strcpy(fullname,cspice_path);
+    strcpy((char*)&(fullname[strlen(cspice_path)]),kname);
+    printf("loading %s\n",fullname);
+    furnsh_c(fullname);
+    free(fullname);
+
+    kname="/moon_de440_220930.tf\0";
+    fullname=(char*)calloc((size_t)strlen((char*)cspice_path)+1+strlen((char*)kname),sizeof(char));
+    if (fullname == 0 ) {
+         fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
+         exit(1);
+    }
+    strcpy(fullname,cspice_path);
+    strcpy((char*)&(fullname[strlen(cspice_path)]),kname);
+    printf("loading %s\n",fullname);
+    furnsh_c(fullname);
+    free(fullname);
+
+    kname="/moon_pa_de440_200625.bpc\0";
+    fullname=(char*)calloc((size_t)strlen((char*)cspice_path)+1+strlen((char*)kname),sizeof(char));
+    if (fullname == 0 ) {
+         fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
+         exit(1);
+    }
+    strcpy(fullname,cspice_path);
+    strcpy((char*)&(fullname[strlen(cspice_path)]),kname);
+    printf("loading %s\n",fullname);
+    furnsh_c(fullname);
+    free(fullname);
+  } else {
+    fprintf(stderr,"CSPICE kernel path 'CSPICE_KERNEL_PATH' is not found in environment variables\n");
+    fprintf(stderr,"Download the kernels pck00011.tpc, naif0012.tls,\n moon_de440_220930.tf, moon_pa_de440_200625.bpc\n");
+    fprintf(stderr,"And rerun the program after setting the directory where these kernels stored as CSPICE_KERNEL_PATH\n");
+    exit(1);
+  }
+}
+
+void
+print_help(void) {
+   fprintf(stderr,"Calculate the UVW coordinates for the given MS based on lunar frame.\n");
+   fprintf(stderr,"Usage:\n");
+   fprintf(stderr,"uvwriter -d MS\n");
+   fprintf(stderr,"-d : input MS (TIME and ANTENNA positions will be used to calculate the UVW coordinates)\n");
+}
+
+
+/* for getopt() */
+extern char *optarg;
+extern int optind, opterr, optopt;
+
+using namespace casacore;
 
 int 
 main(int argc, char **argv) {
 
+  int c;
+  char *inms=0;
+  while ((c=getopt(argc,argv,"d:h"))!=-1) {
+    switch(c) {
+    case 'd':
+      if (optarg) {
+        inms=(char*)calloc((size_t)strlen((char*)optarg)+1,sizeof(char));
+        if (inms== 0 ) {
+         fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
+         exit(1);
+        }
+        strcpy(inms,(char*)optarg);
+      }
+      break;
+    default:
+      print_help();
+      break;
+    }
+  }
+  if (!inms) {
+    print_help();
+    exit(0);
+  }
+
+  cspice_load_kernels();
+
+  MeasurementSet ms(inms,Table::Update);
+  MSAntenna _ant = ms.antenna();
+  size_t N=_ant.nrow();
+  ArrayColumn<double> position(_ant,MSAntenna::columnName(MSAntennaEnums::POSITION));
+  double *xyz=new double[3*N];
+  for (size_t ci=0; ci<N; ci++) {
+    double *_pos=position(ci).data();
+    //printf("Ant %ld %lf %lf %lf\n",ci,_pos[0],_pos[1],_pos[2]);
+    xyz[3*ci]=_pos[0];
+    xyz[3*ci+1]=_pos[1];
+    xyz[3*ci+2]=_pos[2];
+  }
+
+  MSField _field = Table(ms.field());
+  ArrayColumn<double> ref_dir(_field, MSField::columnName(MSFieldEnums::PHASE_DIR));
+  Array<double> dir = ref_dir(0);
+  double *ph = dir.data();
+  double ra0=ph[0];
+  double dec0=ph[1];
+
+  printf("Antennas %ld phase center %lf,%lf (rad)\n",N,ra0,dec0);
+
+  Block<int> sort(1);
+  sort[0]=MS::TIME;
+  MSIter mi(ms,sort,100.0);
+
+  mi.origin();
+  while(mi.more()) {
+    /**************************************************************/
+    Block<String> iv1(3);
+    iv1[0] = "TIME";
+    iv1[1] = "ANTENNA1";
+    iv1[2] = "ANTENNA2";
+    Table t=mi.table().sort(iv1,Sort::Ascending);
+    size_t n_row=t.nrow();
+
+    ROScalarColumn<int> a1(t,"ANTENNA1"), a2(t,"ANTENNA2");
+    ArrayColumn<double> uvwCol(t,"UVW");
+    ROScalarColumn<double> tut(t,"TIME");
+
+    double old_t0=0.0;
+
+    /* unit vector pointing to phase center */
+    double e[3]={0.0,0.0,0.0};
+    for (size_t row=0; row<n_row; row++) {
+      size_t i=a1(row);
+      size_t j=a2(row);
+      Array<double> uvw=uvwCol(row);
+      double *uvwp=uvw.data();
+//      printf("utc %lf ant %ld %ld uvw %lf %lf %lf\n",tut(row),i,j,uvwp[0],uvwp[1],uvwp[2]);
+      double t0=tut(row);
+      if (old_t0 != t0) {
+        /* rotate direction vector to this time */
+        /* convert t0 (s) to JD (days) */
+        double t0_d=t0/86400.0+2400000.5;
+        /* ephemeris time from t0 (s) */
+        double ep_t0=unitim_c(t0_d,"JED","ET");
+        //printf("UTC %le (s) ET %le (s)\n",t0,ep_t0);
+        
+        /* rectangular coords of phace center */
+        SpiceDouble srcrect[3],mtrans[3][3],v2000[3];
+        /* ra,dec to rectangular */
+        radrec_c(1.0,ra0*rpd_c(),dec0*rpd_c(),v2000);
+        /* precess ep_t0 on lunar frame ME: mean Earth/polar axis, PA: principle axis */
+        pxform_c("J2000","MOON_PA",ep_t0,mtrans);//MOON_PA,MOON_ME,IAU_EARTH,ITRF93
+        /* rotate v2000 onto lunar frame */
+        mxv_c(mtrans,v2000,srcrect);
+
+        //printf("source J2000 %lf,%lf,%lf MOON %lf,%lf,%lf\n",v2000[0],v2000[1],v2000[2],
+        // srcrect[0],srcrect[1],srcrect[2]);
+
+        /* unit vector pointing to phase center */
+        double smag=sqrt(srcrect[0]*srcrect[0]+srcrect[1]*srcrect[1]+srcrect[2]*srcrect[2]);
+        e[0]=srcrect[0]/smag;
+        e[1]=srcrect[1]/smag;
+        e[2]=srcrect[2]/smag;
+
+        old_t0=t0;
+      }
+
+      /* find vector stat1-stat2 */
+      if (i!=j) {
+       double v[3];
+       v[0]=xyz[3*i]-xyz[3*j];
+       v[1]=xyz[3*i+1]-xyz[3*j+1];
+       v[2]=xyz[3*i+2]-xyz[3*j+2];
+      
+       /* project v onto plane normal to unit vector e */
+       /* dot product v.e */
+       double v_e=v[0]*e[0]+v[1]*e[1]+v[2]*e[2];
+       /* projected v = v - (v . e) e */
+       v[0] -= v_e*e[0];
+       v[1] -= v_e*e[1];
+       v[2] -= v_e*e[2];
+
+       uvwp[0]=v[0];
+       uvwp[1]=v[1];
+       uvwp[2]=v[2];
+       uvwCol.put(row,uvw);
+      }
+    }
+    /**************************************************************/
+    mi++;
+  }
+
+  delete [] xyz;
+  if (inms) free(inms);
   return 0;
 }

From 208cfc70ec214e305b5e684270c7630078ab5dc6 Mon Sep 17 00:00:00 2001
From: Sarod Yatawatta <sarod@dop373.nfra.nl>
Date: Sun, 21 Apr 2024 09:40:39 +0200
Subject: [PATCH 04/15] working ver

---
 src/uvwriter/uvwriter.cpp | 64 +++++++++++++++++++++++++--------------
 1 file changed, 42 insertions(+), 22 deletions(-)

diff --git a/src/uvwriter/uvwriter.cpp b/src/uvwriter/uvwriter.cpp
index b45cbe5..515a9cb 100644
--- a/src/uvwriter/uvwriter.cpp
+++ b/src/uvwriter/uvwriter.cpp
@@ -91,6 +91,22 @@ cspice_load_kernels(void) {
     printf("loading %s\n",fullname);
     furnsh_c(fullname);
     free(fullname);
+
+    /* following kernel only needed for ITRF93 frame */
+    /*
+    kname="/earth_000101_240713_240419.bpc\0";
+    fullname=(char*)calloc((size_t)strlen((char*)cspice_path)+1+strlen((char*)kname),sizeof(char));
+    if (fullname == 0 ) {
+         fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
+         exit(1);
+    }
+    strcpy(fullname,cspice_path);
+    strcpy((char*)&(fullname[strlen(cspice_path)]),kname);
+    printf("loading %s\n",fullname);
+    furnsh_c(fullname);
+    free(fullname);
+    */
+
   } else {
     fprintf(stderr,"CSPICE kernel path 'CSPICE_KERNEL_PATH' is not found in environment variables\n");
     fprintf(stderr,"Download the kernels pck00011.tpc, naif0012.tls,\n moon_de440_220930.tf, moon_pa_de440_200625.bpc\n");
@@ -149,11 +165,12 @@ main(int argc, char **argv) {
   ArrayColumn<double> position(_ant,MSAntenna::columnName(MSAntennaEnums::POSITION));
   double *xyz=new double[3*N];
   for (size_t ci=0; ci<N; ci++) {
-    double *_pos=position(ci).data();
-    //printf("Ant %ld %lf %lf %lf\n",ci,_pos[0],_pos[1],_pos[2]);
-    xyz[3*ci]=_pos[0];
-    xyz[3*ci+1]=_pos[1];
-    xyz[3*ci+2]=_pos[2];
+    Array<double> _pos=position(ci);
+    double *pos_p=_pos.data();
+    //printf("Ant %ld %lf %lf %lf\n",ci,pos_p[0],pos_p[1],pos_p[2]);
+    xyz[3*ci]=pos_p[0];
+    xyz[3*ci+1]=pos_p[1];
+    xyz[3*ci+2]=pos_p[2];
   }
 
   MSField _field = Table(ms.field());
@@ -185,8 +202,6 @@ main(int argc, char **argv) {
 
     double old_t0=0.0;
 
-    /* unit vector pointing to phase center */
-    double e[3]={0.0,0.0,0.0};
     for (size_t row=0; row<n_row; row++) {
       size_t i=a1(row);
       size_t j=a2(row);
@@ -194,6 +209,7 @@ main(int argc, char **argv) {
       double *uvwp=uvw.data();
 //      printf("utc %lf ant %ld %ld uvw %lf %lf %lf\n",tut(row),i,j,uvwp[0],uvwp[1],uvwp[2]);
       double t0=tut(row);
+      double rotmat[3][3];
       if (old_t0 != t0) {
         /* rotate direction vector to this time */
         /* convert t0 (s) to JD (days) */
@@ -211,15 +227,19 @@ main(int argc, char **argv) {
         /* rotate v2000 onto lunar frame */
         mxv_c(mtrans,v2000,srcrect);
 
+        SpiceDouble s_radius,s_lon,s_lat;
+        reclat_c( srcrect, &s_radius, &s_lon, &s_lat );
+        /* [u,v,w]^T=[sinH cosH 0; -sindel*cosH sindel*sinH cosdel; cosdel*cosH -cosdel*sinH sindel] [x y z]^T */
+        double H=s_lon;
+        double del=M_PI_2-s_lat;
+        rotmat[0][0]=sin(H); rotmat[0][1]=cos(H); rotmat[0][2]=0.0;
+        rotmat[1][0]=-sin(del)*cos(H); rotmat[1][1]=sin(del)*sin(H); rotmat[1][2]=cos(del);
+        rotmat[2][0]=cos(del)*cos(H); rotmat[2][1]=-cos(del)*sin(H); rotmat[2][2]=sin(del);
+
+
         //printf("source J2000 %lf,%lf,%lf MOON %lf,%lf,%lf\n",v2000[0],v2000[1],v2000[2],
         // srcrect[0],srcrect[1],srcrect[2]);
 
-        /* unit vector pointing to phase center */
-        double smag=sqrt(srcrect[0]*srcrect[0]+srcrect[1]*srcrect[1]+srcrect[2]*srcrect[2]);
-        e[0]=srcrect[0]/smag;
-        e[1]=srcrect[1]/smag;
-        e[2]=srcrect[2]/smag;
-
         old_t0=t0;
       }
 
@@ -232,15 +252,15 @@ main(int argc, char **argv) {
       
        /* project v onto plane normal to unit vector e */
        /* dot product v.e */
-       double v_e=v[0]*e[0]+v[1]*e[1]+v[2]*e[2];
-       /* projected v = v - (v . e) e */
-       v[0] -= v_e*e[0];
-       v[1] -= v_e*e[1];
-       v[2] -= v_e*e[2];
-
-       uvwp[0]=v[0];
-       uvwp[1]=v[1];
-       uvwp[2]=v[2];
+       double vloc[3];
+       mxv_c(rotmat,v,vloc);
+
+       //printf("old %lf %lf %lf new %lf %lf %lf\n",uvwp[0],uvwp[1],uvwp[2],vloc[0],vloc[1],vloc[2]);
+
+       uvwp[0]=vloc[0];
+       uvwp[1]=vloc[1];
+       uvwp[2]=vloc[2];
+
        uvwCol.put(row,uvw);
       }
     }

From ad90193ae10c5b4ebe2a319b2f9a991fbf08344a Mon Sep 17 00:00:00 2001
From: Sarod Yatawatta <sarod@dop373.nfra.nl>
Date: Sun, 21 Apr 2024 10:47:38 +0200
Subject: [PATCH 05/15] add cspice lib routines

---
 src/lib/Radio/CMakeLists.txt |   5 ++
 src/lib/Radio/Dirac_radio.h  |   6 ++
 src/lib/Radio/cspice_utils.c | 103 +++++++++++++++++++++++++++++++++++
 src/uvwriter/uvwriter.cpp    |  76 --------------------------
 4 files changed, 114 insertions(+), 76 deletions(-)
 create mode 100644 src/lib/Radio/cspice_utils.c

diff --git a/src/lib/Radio/CMakeLists.txt b/src/lib/Radio/CMakeLists.txt
index 0f393bf..83826e6 100644
--- a/src/lib/Radio/CMakeLists.txt
+++ b/src/lib/Radio/CMakeLists.txt
@@ -14,6 +14,11 @@ set (objects
     diffuse_predict
 )
 
+if(CSPICE_FOUND)
+    set(extra_objects cspice_utils)
+    set(objects ${objects} ${extra_objects})
+endif()
+
 include_directories(${CMAKE_CURRENT_SOURCE_DIR})
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/../Dirac)
 include_directories(${GLIB_PKG_INCLUDE_DIRS})
diff --git a/src/lib/Radio/Dirac_radio.h b/src/lib/Radio/Dirac_radio.h
index 92002a2..0c6db0d 100644
--- a/src/lib/Radio/Dirac_radio.h
+++ b/src/lib/Radio/Dirac_radio.h
@@ -376,6 +376,12 @@ extern int
 sharmonic_modes(int n0,double *th, double *ph, int Nt, complex double *output);
 
 
+/****************************** cspice_utils.c ************************************/
+#ifdef HAVE_CSPICE
+extern void
+cspice_load_kernels(void);
+#endif /* HAVE_CSPICE */
+
 /****************************** shapelet.c ****************************/
 extern complex double
 shapelet_contrib(void*dd, double u, double v, double w);
diff --git a/src/lib/Radio/cspice_utils.c b/src/lib/Radio/cspice_utils.c
new file mode 100644
index 0000000..744fe5b
--- /dev/null
+++ b/src/lib/Radio/cspice_utils.c
@@ -0,0 +1,103 @@
+/*
+ *
+ Copyright (C) 2024 Sarod Yatawatta <sarod@users.sf.net>
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ $Id$
+ */
+
+
+#include <unistd.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include <SpiceUsr.h>
+#include "Dirac_radio.h"
+
+void
+cspice_load_kernels(void) {
+  /* get env CSPICE_KERNEL_PATH */
+  char* cspice_path = getenv("CSPICE_KERNEL_PATH");
+  if (cspice_path) {
+    const char *kname="/pck00011.tpc\0";
+    char *fullname=(char*)calloc((size_t)strlen((char*)cspice_path)+1+strlen((char*)kname),sizeof(char));
+    if (fullname == 0 ) {
+         fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
+         exit(1);
+    }
+    strcpy(fullname,cspice_path);
+    strcpy((char*)&(fullname[strlen(cspice_path)]),kname);
+    printf("loading %s\n",fullname);
+    furnsh_c(fullname);
+    free(fullname);
+
+    kname="/naif0012.tls\0";
+    fullname=(char*)calloc((size_t)strlen((char*)cspice_path)+1+strlen((char*)kname),sizeof(char));
+    if (fullname == 0 ) {
+         fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
+         exit(1);
+    }
+    strcpy(fullname,cspice_path);
+    strcpy((char*)&(fullname[strlen(cspice_path)]),kname);
+    printf("loading %s\n",fullname);
+    furnsh_c(fullname);
+    free(fullname);
+
+    kname="/moon_de440_220930.tf\0";
+    fullname=(char*)calloc((size_t)strlen((char*)cspice_path)+1+strlen((char*)kname),sizeof(char));
+    if (fullname == 0 ) {
+         fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
+         exit(1);
+    }
+    strcpy(fullname,cspice_path);
+    strcpy((char*)&(fullname[strlen(cspice_path)]),kname);
+    printf("loading %s\n",fullname);
+    furnsh_c(fullname);
+    free(fullname);
+
+    kname="/moon_pa_de440_200625.bpc\0";
+    fullname=(char*)calloc((size_t)strlen((char*)cspice_path)+1+strlen((char*)kname),sizeof(char));
+    if (fullname == 0 ) {
+         fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
+         exit(1);
+    }
+    strcpy(fullname,cspice_path);
+    strcpy((char*)&(fullname[strlen(cspice_path)]),kname);
+    printf("loading %s\n",fullname);
+    furnsh_c(fullname);
+    free(fullname);
+
+    /* following kernel only needed for ITRF93 frame */
+    /*
+    kname="/earth_000101_240713_240419.bpc\0";
+    fullname=(char*)calloc((size_t)strlen((char*)cspice_path)+1+strlen((char*)kname),sizeof(char));
+    if (fullname == 0 ) {
+         fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
+         exit(1);
+    }
+    strcpy(fullname,cspice_path);
+    strcpy((char*)&(fullname[strlen(cspice_path)]),kname);
+    printf("loading %s\n",fullname);
+    furnsh_c(fullname);
+    free(fullname);
+    */
+
+  } else {
+    fprintf(stderr,"CSPICE kernel path 'CSPICE_KERNEL_PATH' is not found in environment variables\n");
+    fprintf(stderr,"Download the kernels pck00011.tpc, naif0012.tls,\n moon_de440_220930.tf, moon_pa_de440_200625.bpc\n");
+    fprintf(stderr,"And rerun the program after setting the directory where these kernels stored as CSPICE_KERNEL_PATH\n");
+    exit(1);
+  }
+}
diff --git a/src/uvwriter/uvwriter.cpp b/src/uvwriter/uvwriter.cpp
index 515a9cb..aa53ae9 100644
--- a/src/uvwriter/uvwriter.cpp
+++ b/src/uvwriter/uvwriter.cpp
@@ -39,82 +39,6 @@
 
 //#define DEBUG
 
-void
-cspice_load_kernels(void) {
-  /* get env CSPICE_KERNEL_PATH */
-  char* cspice_path = getenv("CSPICE_KERNEL_PATH");
-  if (cspice_path) {
-    const char *kname="/pck00011.tpc\0";
-    char *fullname=(char*)calloc((size_t)strlen((char*)cspice_path)+1+strlen((char*)kname),sizeof(char));
-    if (fullname == 0 ) {
-         fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
-         exit(1);
-    }
-    strcpy(fullname,cspice_path);
-    strcpy((char*)&(fullname[strlen(cspice_path)]),kname);
-    printf("loading %s\n",fullname);
-    furnsh_c(fullname);
-    free(fullname);
-
-    kname="/naif0012.tls\0";
-    fullname=(char*)calloc((size_t)strlen((char*)cspice_path)+1+strlen((char*)kname),sizeof(char));
-    if (fullname == 0 ) {
-         fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
-         exit(1);
-    }
-    strcpy(fullname,cspice_path);
-    strcpy((char*)&(fullname[strlen(cspice_path)]),kname);
-    printf("loading %s\n",fullname);
-    furnsh_c(fullname);
-    free(fullname);
-
-    kname="/moon_de440_220930.tf\0";
-    fullname=(char*)calloc((size_t)strlen((char*)cspice_path)+1+strlen((char*)kname),sizeof(char));
-    if (fullname == 0 ) {
-         fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
-         exit(1);
-    }
-    strcpy(fullname,cspice_path);
-    strcpy((char*)&(fullname[strlen(cspice_path)]),kname);
-    printf("loading %s\n",fullname);
-    furnsh_c(fullname);
-    free(fullname);
-
-    kname="/moon_pa_de440_200625.bpc\0";
-    fullname=(char*)calloc((size_t)strlen((char*)cspice_path)+1+strlen((char*)kname),sizeof(char));
-    if (fullname == 0 ) {
-         fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
-         exit(1);
-    }
-    strcpy(fullname,cspice_path);
-    strcpy((char*)&(fullname[strlen(cspice_path)]),kname);
-    printf("loading %s\n",fullname);
-    furnsh_c(fullname);
-    free(fullname);
-
-    /* following kernel only needed for ITRF93 frame */
-    /*
-    kname="/earth_000101_240713_240419.bpc\0";
-    fullname=(char*)calloc((size_t)strlen((char*)cspice_path)+1+strlen((char*)kname),sizeof(char));
-    if (fullname == 0 ) {
-         fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
-         exit(1);
-    }
-    strcpy(fullname,cspice_path);
-    strcpy((char*)&(fullname[strlen(cspice_path)]),kname);
-    printf("loading %s\n",fullname);
-    furnsh_c(fullname);
-    free(fullname);
-    */
-
-  } else {
-    fprintf(stderr,"CSPICE kernel path 'CSPICE_KERNEL_PATH' is not found in environment variables\n");
-    fprintf(stderr,"Download the kernels pck00011.tpc, naif0012.tls,\n moon_de440_220930.tf, moon_pa_de440_200625.bpc\n");
-    fprintf(stderr,"And rerun the program after setting the directory where these kernels stored as CSPICE_KERNEL_PATH\n");
-    exit(1);
-  }
-}
-
 void
 print_help(void) {
    fprintf(stderr,"Calculate the UVW coordinates for the given MS based on lunar frame.\n");

From 97110781b005b90c6e79c5a40a6c59e35c669c88 Mon Sep 17 00:00:00 2001
From: Sarod Yatawatta <sarod@dop373.nfra.nl>
Date: Sun, 21 Apr 2024 10:53:52 +0200
Subject: [PATCH 06/15] fix build

---
 src/MPI/CMakeLists.txt | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/MPI/CMakeLists.txt b/src/MPI/CMakeLists.txt
index 94a92b1..6d2c1cd 100644
--- a/src/MPI/CMakeLists.txt
+++ b/src/MPI/CMakeLists.txt
@@ -41,6 +41,7 @@ if(HAVE_CUDA)
     ${LAPACK_LIBRARIES}
     ${WCSLIB_LIBRARIES}
     ${GLIB_PKG_LIBRARIES}
+    $<$<1:${CSPICE_FOUND}>:${CSPICE_LINK_LIBRARIES}>
     ${MPI_CXX_LIBRARIES}
     ${MPI_CXX_LINK_FLAGS}
     -lpthread
@@ -66,6 +67,7 @@ else()
     ${LAPACK_LIBRARIES}
     ${WCSLIB_LIBRARIES}
     ${GLIB_PKG_LIBRARIES}
+    $<$<1:${CSPICE_FOUND}>:${CSPICE_LINK_LIBRARIES}>
     ${MPI_CXX_LIBRARIES}
     ${MPI_CXX_LINK_FLAGS}
     -lpthread

From b8d2c79fd189276ac55c99b48608ce6de57d5091 Mon Sep 17 00:00:00 2001
From: Sarod Yatawatta <sarod@dop373.nfra.nl>
Date: Sun, 21 Apr 2024 11:29:26 +0200
Subject: [PATCH 07/15] fix build

---
 src/MPI/CMakeLists.txt      | 4 ++--
 src/uvwriter/CMakeLists.txt | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/MPI/CMakeLists.txt b/src/MPI/CMakeLists.txt
index 6d2c1cd..9a0a185 100644
--- a/src/MPI/CMakeLists.txt
+++ b/src/MPI/CMakeLists.txt
@@ -41,7 +41,7 @@ if(HAVE_CUDA)
     ${LAPACK_LIBRARIES}
     ${WCSLIB_LIBRARIES}
     ${GLIB_PKG_LIBRARIES}
-    $<$<1:${CSPICE_FOUND}>:${CSPICE_LINK_LIBRARIES}>
+    $<$<BOOL:${CSPICE_FOUND}>:${CSPICE_LINK_LIBRARIES}>
     ${MPI_CXX_LIBRARIES}
     ${MPI_CXX_LINK_FLAGS}
     -lpthread
@@ -67,7 +67,7 @@ else()
     ${LAPACK_LIBRARIES}
     ${WCSLIB_LIBRARIES}
     ${GLIB_PKG_LIBRARIES}
-    $<$<1:${CSPICE_FOUND}>:${CSPICE_LINK_LIBRARIES}>
+    $<$<BOOL:${CSPICE_FOUND}>:${CSPICE_LINK_LIBRARIES}>
     ${MPI_CXX_LIBRARIES}
     ${MPI_CXX_LINK_FLAGS}
     -lpthread
diff --git a/src/uvwriter/CMakeLists.txt b/src/uvwriter/CMakeLists.txt
index 5949c03..e0ceaa5 100644
--- a/src/uvwriter/CMakeLists.txt
+++ b/src/uvwriter/CMakeLists.txt
@@ -39,7 +39,7 @@ if (HAVE_CUDA)
         ${LAPACK_LIBRARIES}
         ${WCSLIB_LIBRARIES}
         ${GLIB_PKG_LIBRARIES}
-        $<$<1:${CSPICE_FOUND}>:${CSPICE_LINK_LIBRARIES}>
+        $<$<BOOL:${CSPICE_FOUND}>:${CSPICE_LINK_LIBRARIES}>
         -lpthread
         -lm
         ${CUDA_CUBLAS_LIBRARIES}
@@ -60,7 +60,7 @@ else()
         ${LAPACK_LIBRARIES}
         ${WCSLIB_LIBRARIES}
         ${GLIB_PKG_LIBRARIES}
-        $<$<1:${CSPICE_FOUND}>:${CSPICE_LINK_LIBRARIES}>
+        $<$<BOOL:${CSPICE_FOUND}>:${CSPICE_LINK_LIBRARIES}>
         -lpthread
         -lm
   )

From 9a0f9dac13240a2258f2338a7db6f118d065613e Mon Sep 17 00:00:00 2001
From: Sarod Yatawatta <sarod@dop373.nfra.nl>
Date: Sun, 21 Apr 2024 11:48:04 +0200
Subject: [PATCH 08/15] add links

---
 src/lib/Radio/cspice_utils.c | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/lib/Radio/cspice_utils.c b/src/lib/Radio/cspice_utils.c
index 744fe5b..119016d 100644
--- a/src/lib/Radio/cspice_utils.c
+++ b/src/lib/Radio/cspice_utils.c
@@ -96,7 +96,11 @@ cspice_load_kernels(void) {
 
   } else {
     fprintf(stderr,"CSPICE kernel path 'CSPICE_KERNEL_PATH' is not found in environment variables\n");
-    fprintf(stderr,"Download the kernels pck00011.tpc, naif0012.tls,\n moon_de440_220930.tf, moon_pa_de440_200625.bpc\n");
+    fprintf(stderr,"Download the kernels\n"
+    "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/pck00011.tpc\n"
+    "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls\n"
+    "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/fk/satellites/moon_de440_220930.tf\n"
+    "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/moon_pa_de440_200625.bpc\n");
     fprintf(stderr,"And rerun the program after setting the directory where these kernels stored as CSPICE_KERNEL_PATH\n");
     exit(1);
   }

From 30e6faf1fca800e68fd143a41c1447f992d77bd5 Mon Sep 17 00:00:00 2001
From: Sarod Yatawatta <sarod@dop373.nfra.nl>
Date: Tue, 23 Apr 2024 12:06:42 +0200
Subject: [PATCH 09/15] in lunar mode no precession

---
 src/MPI/sagecal_slave.cpp            | 2 +-
 src/MPI/sagecal_stochastic_slave.cpp | 2 +-
 src/MS/fullbatch_mode.cpp            | 2 +-
 src/MS/minibatch_consensus_mode.cpp  | 2 +-
 src/MS/minibatch_mode.cpp            | 2 +-
 5 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/MPI/sagecal_slave.cpp b/src/MPI/sagecal_slave.cpp
index 5bb1939..e494723 100644
--- a/src/MPI/sagecal_slave.cpp
+++ b/src/MPI/sagecal_slave.cpp
@@ -618,7 +618,7 @@ cout<<myrank<<" : "<<cm<<": downweight ratio ("<<iodata_vec[cm].fratio<<") based
      for(int cm=0; cm<mymscount; cm++) {
      /* precess source locations (also beam pointing) from J2000 to JAPP if we do any beam predictions,
       using first time slot as epoch */
-      if (doBeam && !sources_precessed) {
+      if (doBeam && !sources_precessed && (beam_vec[cm].elType!=ELEM_ALO)) {
         Data::precess_source_locations(beam_vec[cm].time_utc[iodata_vec[cm].tilesz/2],carr_vec[cm],M,&beam_vec[cm].p_ra0,&beam_vec[cm].p_dec0,&beam_vec[cm].b_ra0,&beam_vec[cm].b_dec0,Data::Nt);
        sources_precessed=1;
       }
diff --git a/src/MPI/sagecal_stochastic_slave.cpp b/src/MPI/sagecal_stochastic_slave.cpp
index 0c1d4ea..f7a86d3 100644
--- a/src/MPI/sagecal_stochastic_slave.cpp
+++ b/src/MPI/sagecal_stochastic_slave.cpp
@@ -692,7 +692,7 @@ cout<<"Worker "<<myrank<<" quitting"<<endl;
 
         /* precess source locations (also beam pointing) from J2000 to JAPP if we do any beam predictions,
            using first time slot as epoch */
-        if (doBeam && !sources_precessed) {
+        if (doBeam && !sources_precessed && (beam_vec[cm].elType!=ELEM_ALO)) {
           Data::precess_source_locations(beam_vec[cm].time_utc[iodata_vec[cm].tilesz/2],carr_vec[cm],M,&beam_vec[cm].p_ra0,&beam_vec[cm].p_dec0,&beam_vec[cm].b_ra0,&beam_vec[cm].b_dec0,Data::Nt);
          if (cm==mymscount-1) {sources_precessed=1;} /* wait till all MS are handled */
         }
diff --git a/src/MS/fullbatch_mode.cpp b/src/MS/fullbatch_mode.cpp
index d2b3a36..a1a7083 100644
--- a/src/MS/fullbatch_mode.cpp
+++ b/src/MS/fullbatch_mode.cpp
@@ -324,7 +324,7 @@ run_fullbatch_calibration(void) {
     }
     /* precess source locations (also beam pointing) from J2000 to JAPP if we do any beam predictions,
       using first time slot as epoch */
-    if (doBeam && !sources_precessed) {
+    if (doBeam && !sources_precessed && (beam.elType!=ELEM_ALO)) {
       Data::precess_source_locations(beam.time_utc[iodata.tilesz/2],carr,M,&beam.p_ra0,&beam.p_dec0,&beam.b_ra0,&beam.b_dec0,Data::Nt);
       sources_precessed=1;
     }
diff --git a/src/MS/minibatch_consensus_mode.cpp b/src/MS/minibatch_consensus_mode.cpp
index bc89106..167d920 100644
--- a/src/MS/minibatch_consensus_mode.cpp
+++ b/src/MS/minibatch_consensus_mode.cpp
@@ -479,7 +479,7 @@ run_minibatch_consensus_calibration(void) {
 
         /* precess source locations (also beam pointing) from J2000 to JAPP if we do any beam predictions,
            using first time slot as epoch */
-        if (doBeam && !sources_precessed) {
+        if (doBeam && !sources_precessed && (beam.elType!=ELEM_ALO)) {
          Data::precess_source_locations(beam.time_utc[iodata.tilesz/2],carr,M,&beam.p_ra0,&beam.p_dec0,&beam.b_ra0,&beam.b_dec0,Data::Nt);
          sources_precessed=1;
         }
diff --git a/src/MS/minibatch_mode.cpp b/src/MS/minibatch_mode.cpp
index f2bc5cc..ee1414b 100644
--- a/src/MS/minibatch_mode.cpp
+++ b/src/MS/minibatch_mode.cpp
@@ -391,7 +391,7 @@ run_minibatch_calibration(void) {
 
         /* precess source locations (also beam pointing) from J2000 to JAPP if we do any beam predictions,
            using first time slot as epoch */
-        if (doBeam && !sources_precessed) {
+        if (doBeam && !sources_precessed && (beam.elType!=ELEM_ALO)) {
           Data::precess_source_locations(beam.time_utc[iodata.tilesz/2],carr,M,&beam.p_ra0,&beam.p_dec0,&beam.b_ra0,&beam.b_dec0,Data::Nt);
          sources_precessed=1;
         }

From 46ef31513597b9abdebc9d1c258ac3803438ae8b Mon Sep 17 00:00:00 2001
From: Sarod Yatawatta <sarod@dop373.nfra.nl>
Date: Fri, 26 Apr 2024 23:13:20 +0200
Subject: [PATCH 10/15] add cspice build scripts

---
 scripts/CSPICE/README.md                 |  39 +++++++
 scripts/CSPICE/makeall.sh                | 139 +++++++++++++++++++++++
 scripts/CSPICE/mkprodct.sh               | 132 +++++++++++++++++++++
 scripts/CSPICE/pkgconfig.stub            |  12 ++
 scripts/load_das5_astron_modules_gcc4.sh |  14 ---
 scripts/load_das5_astron_modules_gcc6.sh |  14 ---
 scripts/sagecal_gpu.job                  |  26 -----
 7 files changed, 322 insertions(+), 54 deletions(-)
 create mode 100644 scripts/CSPICE/README.md
 create mode 100755 scripts/CSPICE/makeall.sh
 create mode 100755 scripts/CSPICE/mkprodct.sh
 create mode 100644 scripts/CSPICE/pkgconfig.stub
 delete mode 100755 scripts/load_das5_astron_modules_gcc4.sh
 delete mode 100755 scripts/load_das5_astron_modules_gcc6.sh
 delete mode 100644 scripts/sagecal_gpu.job

diff --git a/scripts/CSPICE/README.md b/scripts/CSPICE/README.md
new file mode 100644
index 0000000..22e07a5
--- /dev/null
+++ b/scripts/CSPICE/README.md
@@ -0,0 +1,39 @@
+# Installing CSPICE toolkit
+This document describes building the CSPICE toolkit to link with sagecal.
+
+# Download
+Get the latest [cspice.tar.Z](https://naif.jpl.nasa.gov/pub/naif/toolkit//C/PC_Linux_GCC_64bit/packages/cspice.tar.Z)
+
+```
+wget https://naif.jpl.nasa.gov/pub/naif/toolkit//C/PC_Linux_GCC_64bit/packages/cspice.tar.Z
+```
+
+Extract the files
+
+```
+tar xvf cspice.tar.Z
+```
+and there will be directory ```cspice```.
+
+# Copy files
+Copy all files in this directory to the ```cspice``` directory.
+
+```
+cp makeall.sh mkprodct.sh pkgconfig.stub /full_path_to/cspice/
+```
+
+Note that ```/full_path_to/cspice/``` is the directory created by extracting cspice.tar.Z.
+
+# Build
+In the ```cspice``` directory, run
+
+```
+bash ./makeall.sh
+```
+and you are ready to link with sagecal.
+
+
+# Linking with sagecal
+Use cmake flag ```-DCSPICE_PREFIX=/full_path_to/cspice/lib``` to point to ```cspice/lib``` directory created after builiding CSPICE.
+
+vr 26 apr 2024 23:08:20 CEST
diff --git a/scripts/CSPICE/makeall.sh b/scripts/CSPICE/makeall.sh
new file mode 100755
index 0000000..a20cb53
--- /dev/null
+++ b/scripts/CSPICE/makeall.sh
@@ -0,0 +1,139 @@
+echo This script builds the SPICE delivery
+echo for the cspice package of the toolkit.
+echo
+echo The script must be executed from the
+echo cspice directory.
+echo
+cd src
+echo
+echo Creating cspice
+echo
+cd cspice
+../../mkprodct.sh
+cd ..
+echo
+echo Creating csupport
+echo
+cd csupport
+../../mkprodct.sh
+cd ..
+echo
+echo Creating brief_c
+echo
+cd brief_c
+../../mkprodct.sh
+cd ..
+echo
+echo Creating chrnos_c
+echo
+cd chrnos_c
+../../mkprodct.sh
+cd ..
+echo
+echo Creating ckbref_c
+echo
+cd ckbref_c
+../../mkprodct.sh
+cd ..
+echo
+echo Creating commnt_c
+echo
+cd commnt_c
+../../mkprodct.sh
+cd ..
+echo
+echo Creating cook_c
+echo
+cd cook_c
+../../mkprodct.sh
+cd ..
+echo
+echo Creating dskbrief_c
+echo
+cd dskbrief_c
+../../mkprodct.sh
+cd ..
+echo
+echo Creating dskexp_c
+echo
+cd dskexp_c
+../../mkprodct.sh
+cd ..
+echo
+echo Creating frmdif_c
+echo
+cd frmdif_c
+../../mkprodct.sh
+cd ..
+echo
+echo Creating inspkt_c
+echo
+cd inspkt_c
+../../mkprodct.sh
+cd ..
+echo
+echo Creating mkdsk_c
+echo
+cd mkdsk_c
+../../mkprodct.sh
+cd ..
+echo
+echo Creating mkspk_c
+echo
+cd mkspk_c
+../../mkprodct.sh
+cd ..
+echo
+echo Creating msopck_c
+echo
+cd msopck_c
+../../mkprodct.sh
+cd ..
+echo
+echo Creating spacit_c
+echo
+cd spacit_c
+../../mkprodct.sh
+cd ..
+echo
+echo Creating spkdif_c
+echo
+cd spkdif_c
+../../mkprodct.sh
+cd ..
+echo
+echo Creating spkmrg_c
+echo
+cd spkmrg_c
+../../mkprodct.sh
+cd ..
+echo
+echo Creating tobin_c
+echo
+cd tobin_c
+../../mkprodct.sh
+cd ..
+echo
+echo Creating toxfr_c
+echo
+cd toxfr_c
+../../mkprodct.sh
+cd ..
+echo
+echo Creating versn_c
+echo
+cd versn_c
+../../mkprodct.sh
+cd ..
+cd ..
+echo Toolkit Build Complete
+
+echo Making links
+cd ./lib && ln -s cspice.a libcspice.a && ln -s csupport.a libcsupport.a && cd ..
+echo Done making links
+
+echo Making pkgconfig file
+mkdir -p ./lib/pkgconfig && echo "prefix=\"`pwd`\"" > ./lib/pkgconfig/libcspice.pc && cat pkgconfig.stub >> ./lib/pkgconfig/libcspice.pc
+echo Done making pkgconfig
+
+echo All done
diff --git a/scripts/CSPICE/mkprodct.sh b/scripts/CSPICE/mkprodct.sh
new file mode 100755
index 0000000..efbd83a
--- /dev/null
+++ b/scripts/CSPICE/mkprodct.sh
@@ -0,0 +1,132 @@
+#!/bin/bash
+
+
+# find .pgm files (if main.x is given)
+if [ -f main.x ]; then
+if compgen -G "*.pgm" > /dev/null; then
+ for MAIN in *.pgm; do
+  echo $MAIN;
+  STEM=`echo $MAIN|sed 's/.pgm//g'`
+  TARGET=$STEM.px
+  cp $MAIN $STEM"_main.c"
+  cp main.x $TARGET
+ done
+fi
+fi
+
+CC=gcc
+if ! [ -f main.x ]; then
+ CFLAGS="-c -ansi -m64 -O2 -DNON_UNIX_STDIO"
+else
+ CFLAGS="-c -ansi -m64 -O2 -fPIC -DNON_UNIX_STDIO"
+fi
+LDFLAGS="-lm -m64"
+# set compiler to gcc
+TKCOMPILER=$CC
+TKCOMPILEOPTIONS=$CFLAGS
+TKLINKOPTIONS=$LDFLAGS
+
+
+# library name is current dir
+DIR=$(basename `pwd`)
+LIBRARY="../../lib/"$DIR
+
+
+# check any *.c files need compiling
+if compgen -G "*.c" > /dev/null; then
+ for SRCFILE in *.c; do
+   $TKCOMPILER $TKCOMPILEOPTIONS $SRCFILE
+ done
+fi
+
+# if program files exit, need to create a library (override LIBRARY above)
+if compgen -G "*.pgm" > /dev/null; then
+  LIBRARY="locallib"
+fi
+
+if compgen -G "*.o" > /dev/null; then
+  echo "linking library $LIBRARY"
+  ar crv $LIBRARY.a *.o
+  ranlib $LIBRARY.a
+  rm *.o
+fi
+
+
+# if there are programs, compile them
+if compgen -G "*.pgm" > /dev/null; then
+  # if main.x is not present, go through *.pgm, not *.px files
+  if ! [ -f main.x ]; then
+    for MAIN in *.pgm; do
+      STEM=`echo $MAIN|sed 's/.pgm//g'`
+      TARGET=$STEM.c
+      MAINOBJ=$STEM.o
+      EXECUT="../../exe/"$STEM
+      cp $MAIN $TARGET
+      echo "compiling and linking $MAIN"
+      if [ -f locallib.a ]; then
+        $TKCOMPILER $TKCOMPILEOPTIONS $TARGET
+        $TKCOMPILER -o $EXECUT $MAINOBJ \
+          locallib.a \
+          ../../lib/csupport.a \
+          ../../lib/cspice.a \
+          $TKLINKOPTIONS
+
+        rm $TARGET
+        rm $MAINOBJ
+        rm locallib.a
+      else
+        $TKCOMPILER $TKCOMPILEOPTIONS $TARGET
+        $TKCOMPILER -o $EXECUT $MAINOBJ \
+          ../../lib/csupport.a \
+          ../../lib/cspice.a \
+          $TKLINKOPTIONS
+
+        rm $TARGET
+        rm $MAINOBJ
+      fi
+    done
+  else
+    for MAIN in *.px; do
+      STEM=`echo $MAIN|sed 's/.px//g'`
+      TARGET=$STEM.c
+      MAINOBJ=$STEM.o
+      EXECUT="../../exe/"$STEM
+      cp $MAIN $TARGET
+      echo "compiling and linking $MAIN"
+      if [ -f locallib.a ]; then
+        $TKCOMPILER $TKCOMPILEOPTIONS $TARGET
+        $TKCOMPILER -o $EXECUT $MAINOBJ \
+          locallib.a \
+          ../../lib/csupport.a \
+          ../../lib/cspice.a \
+          $TKLINKOPTIONS
+
+        rm $TARGET
+        rm $MAINOBJ
+        rm locallib.a
+      else
+        $TKCOMPILER $TKCOMPILEOPTIONS $TARGET
+        $TKCOMPILER -o $EXECUT $MAINOBJ \
+          ../../lib/csupport.a \
+          ../../lib/cspice.a \
+          $TKLINKOPTIONS
+
+        rm $TARGET
+        rm $MAINOBJ
+      fi
+    done
+  fi
+fi
+
+
+
+# cleanup
+if compgen -G "*.o" > /dev/null; then
+ rm *.o
+fi
+if compgen -G "*.px" > /dev/null; then
+ rm *.px
+fi
+if compgen -G "*_main.c" > /dev/null; then
+ rm *_main.c
+fi
diff --git a/scripts/CSPICE/pkgconfig.stub b/scripts/CSPICE/pkgconfig.stub
new file mode 100644
index 0000000..67c5fa2
--- /dev/null
+++ b/scripts/CSPICE/pkgconfig.stub
@@ -0,0 +1,12 @@
+exec_prefix="${prefix}"
+libdir="${prefix}/lib"
+includedir="${prefix}/include"
+
+Name: CSPICE
+Description: 
+URL:
+Version: 0.0.1
+Requires: 
+Requires.private: 
+Cflags: -I"${includedir}"
+Libs: -L"${libdir}" -lcspice
diff --git a/scripts/load_das5_astron_modules_gcc4.sh b/scripts/load_das5_astron_modules_gcc4.sh
deleted file mode 100755
index d4a5f9b..0000000
--- a/scripts/load_das5_astron_modules_gcc4.sh
+++ /dev/null
@@ -1,14 +0,0 @@
-#!/bin/bash
-
-module load cmake/3.8.2
-module load mpich/ge/gcc/64/3.2
-module load gcc/4.9.3
-module load casacore/2.3.0-gcc-4.9.3
-module load wcslib/5.16-gcc-4.9.3
-module load cfitsio/3.410-gcc-4.9.3
-#module load blas/gcc/64/3.7.0
-module load openblas/0.2.20
-
-
-module load cuda91/toolkit/9.1.85
-#module load cuda92/toolkit/9.2.88
diff --git a/scripts/load_das5_astron_modules_gcc6.sh b/scripts/load_das5_astron_modules_gcc6.sh
deleted file mode 100755
index 16dfff7..0000000
--- a/scripts/load_das5_astron_modules_gcc6.sh
+++ /dev/null
@@ -1,14 +0,0 @@
-#!/bin/bash
-
-module load cmake/3.8.2
-module load mpich/ge/gcc/64/3.2
-module load gcc/6.3.0
-module load casacore/2.4.1-gcc-6.3.0
-module load wcslib/5.18-gcc-6.3.0
-module load cfitsio/3.430-gcc-6.3.0
-#module load blas/gcc/64/3.7.0
-module load openblas/0.2.20
-
-
-module load cuda91/toolkit/9.1.85
-#module load cuda92/toolkit/9.2.88
diff --git a/scripts/sagecal_gpu.job b/scripts/sagecal_gpu.job
deleted file mode 100644
index f07d147..0000000
--- a/scripts/sagecal_gpu.job
+++ /dev/null
@@ -1,26 +0,0 @@
-#!/bin/sh
-#SBATCH --time=00:2:00
-#SBATCH -N 1
-###SBATCH -C TitanX
-###SBATCH --gres=gpu:1
-
-. /etc/bashrc
-. /etc/profile.d/modules.sh
-
-module load cmake/3.8.2
-module load mpich/ge/gcc/64/3.2
-module load gcc/6.3.0
-module load casacore/2.4.1-gcc-6.3.0
-module load wcslib/5.18-gcc-6.3.0
-module load cfitsio/3.430-gcc-6.3.0
-#module load blas/gcc/64/3.7.0
-module load openblas/0.2.20
-module load cuda91/toolkit/9.1.85
-
-nvidia-smi
-
-export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/cm/shared/package/cuda91/toolkit/9.1.85/lib64/stubs
-ln -s /cm/shared/package/cuda91/toolkit/9.1.85/lib64/stubs/libnvidia-ml.so libnvidia-ml.so.1
-
-/home/fdiblen/sagecal_gpu/build/dist/bin/sagecal_gpu
-

From cb66fade758e87ebe00aad3cd6d68c58a2b80e80 Mon Sep 17 00:00:00 2001
From: Sarod Yatawatta <sarod@dop373.nfra.nl>
Date: Sun, 28 Apr 2024 18:11:17 +0200
Subject: [PATCH 11/15] full cspice integration

---
 scripts/CSPICE/README.md         |   4 +-
 src/MS/data.cpp                  |  42 +++++++++-
 src/MS/fullbatch_mode.cpp        |  19 +++--
 src/MS/main.cpp                  |   7 ++
 src/lib/Dirac/Dirac_common.h     |  12 +++
 src/lib/Radio/Dirac_radio.h      |   9 ++
 src/lib/Radio/cspice_utils.c     | 139 ++++++++++++++++++++++++++++++-
 src/lib/Radio/predict_withbeam.c |  47 +++++++++--
 8 files changed, 261 insertions(+), 18 deletions(-)

diff --git a/scripts/CSPICE/README.md b/scripts/CSPICE/README.md
index 22e07a5..7774dde 100644
--- a/scripts/CSPICE/README.md
+++ b/scripts/CSPICE/README.md
@@ -34,6 +34,6 @@ and you are ready to link with sagecal.
 
 
 # Linking with sagecal
-Use cmake flag ```-DCSPICE_PREFIX=/full_path_to/cspice/lib``` to point to ```cspice/lib``` directory created after builiding CSPICE.
+Use cmake flag ```-DCSPICE_PREFIX=/full_path_to/cspice``` to point to ```cspice``` directory created after exctracting CSPICE.
 
-vr 26 apr 2024 23:08:20 CEST
+za 27 apr 2024  9:23:45 CEST
diff --git a/src/MS/data.cpp b/src/MS/data.cpp
index e316fc5..261f590 100644
--- a/src/MS/data.cpp
+++ b/src/MS/data.cpp
@@ -236,6 +236,11 @@ Data::readAuxData(const char *fname, Data::IOData *data, Data::LBeam *binfo) {
      binfo->elType=(data->freq0<100e6?ELEM_LBA:ELEM_HBA);
    } else if ( !tel.compare("ALO") ) {
      binfo->elType=ELEM_ALO;
+#ifdef HAVE_CSPICE
+     cspice_load_kernels();
+#else
+     std::cout<<"Warning: telecope "<<tel<<", but CSPICE is not found"<<std::endl;
+#endif
    } else {
      std::cout<<"Warning: unknown telecope "<<tel<<", defaulting to LOFAR"<<std::endl;
      binfo->elType=(data->freq0<100e6?ELEM_LBA:ELEM_HBA);
@@ -291,12 +296,30 @@ Data::readAuxData(const char *fname, Data::IOData *data, Data::LBeam *binfo) {
        double *tx=_pos.data();
        binfo->sz[ci]=tx[2];
 
+       /* for ALO, use CSPICE to calculate longitude, latitude */
+#ifdef HAVE_CSPICE
+       if (binfo->elType==ELEM_ALO) {
+         double lat,lon,alt;
+         cspice_xyz_to_latlon(tx[0],tx[1],tx[2],&lon,&lat,&alt);
+         binfo->sx[ci]=lon;
+         binfo->sy[ci]=lat;
+         binfo->sz[ci]=alt;
+       } else {
+         MPosition stnpos(MVPosition(tx[0],tx[1],tx[2]),MPosition::ITRF);
+         Array<double> _radpos=stnpos.getAngle("rad").getValue();
+         tx=_radpos.data();
+         binfo->sx[ci]=tx[0];
+         binfo->sy[ci]=tx[1];
+       }
+#else
        MPosition stnpos(MVPosition(tx[0],tx[1],tx[2]),MPosition::ITRF);
        Array<double> _radpos=stnpos.getAngle("rad").getValue();
        tx=_radpos.data();
 
        binfo->sx[ci]=tx[0];
        binfo->sy[ci]=tx[1];
+#endif /* HAVE_CSPICE */
+
        /* following is the number of tiles */
        binfo->Nelem[ci]=offset.shape(ci)[1];
      }
@@ -455,12 +478,29 @@ Data::readAuxData(const char *fname, Data::IOData *data, Data::LBeam *binfo) {
        double *tx=_pos.data();
        binfo->sz[ci]=tx[2];
 
-       MPosition stnpos(MVPosition(tx[0],tx[1],tx[2]), MPosition::ITRF);
+       /* for ALO, use CSPICE to calculate longitude, latitude */
+#ifdef HAVE_CSPICE
+       if (binfo->elType==ELEM_ALO) {
+         double lat,lon,alt;
+         cspice_xyz_to_latlon(tx[0],tx[1],tx[2],&lon,&lat,&alt);
+         binfo->sx[ci]=lon;
+         binfo->sy[ci]=lat;
+         binfo->sz[ci]=alt;
+       } else {
+         MPosition stnpos(MVPosition(tx[0],tx[1],tx[2]),MPosition::ITRF);
+         Array<double> _radpos=stnpos.getAngle("rad").getValue();
+         tx=_radpos.data();
+         binfo->sx[ci]=tx[0];
+         binfo->sy[ci]=tx[1];
+       }
+#else
+       MPosition stnpos(MVPosition(tx[0],tx[1],tx[2]),MPosition::ITRF);
        Array<double> _radpos=stnpos.getAngle("rad").getValue();
        tx=_radpos.data();
 
        binfo->sx[ci]=tx[0];
        binfo->sy[ci]=tx[1];
+#endif /* HAVE_CSPICE */
 
        /* allocate storage for only 1 element */
        binfo->Nelem[ci]=1;
diff --git a/src/MS/fullbatch_mode.cpp b/src/MS/fullbatch_mode.cpp
index a1a7083..9307fd9 100644
--- a/src/MS/fullbatch_mode.cpp
+++ b/src/MS/fullbatch_mode.cpp
@@ -67,9 +67,9 @@ run_fullbatch_calibration(void) {
     if (Data::randomize) {
      srand(time(0)); /* use different seed */
     }
-    if (doBeam==DOBEAM_FULL||doBeam==DOBEAM_ELEMENT) {
+    if (doBeam==DOBEAM_FULL||doBeam==DOBEAM_ELEMENT||doBeam==DOBEAM_ALO) {
      set_elementcoeffs(beam.elType, iodata.freq0, &ecoeff);
-    } else if (doBeam==DOBEAM_FULL_WB||doBeam==DOBEAM_ELEMENT_WB) {
+    } else if (doBeam==DOBEAM_FULL_WB||doBeam==DOBEAM_ELEMENT_WB||doBeam==DOBEAM_ALO_WB) {
      set_elementcoeffs_wb(beam.elType, iodata.freqs, iodata.Nchan, &ecoeff);
     }
 
@@ -532,11 +532,15 @@ run_fullbatch_calibration(void) {
       predict_visibilities_multifreq(iodata.u,iodata.v,iodata.w,iodata.xo,iodata.N,iodata.Nbase,iodata.tilesz,barr,carr,M,iodata.freqs,iodata.Nchan,iodata.deltaf,iodata.deltat,iodata.dec0,Data::Nt,Data::DoSim);
      } else {
       predict_visibilities_multifreq_withbeam(iodata.u,iodata.v,iodata.w,iodata.xo,iodata.N,iodata.Nbase,iodata.tilesz,barr,carr,M,iodata.freqs,iodata.Nchan,iodata.deltaf,iodata.deltat,iodata.dec0,
-  beam.bfType,beam.b_ra0,beam.b_dec0,beam.p_ra0,beam.p_dec0,iodata.freq0,beam.sx,beam.sy,beam.time_utc,beam.Nelem,beam.xx,beam.yy,beam.zz,&ecoeff,doBeam,Data::Nt,Data::DoSim);
+      beam.bfType,beam.b_ra0,beam.b_dec0,beam.p_ra0,beam.p_dec0,iodata.freq0,beam.sx,beam.sy,beam.time_utc,beam.Nelem,beam.xx,beam.yy,beam.zz,&ecoeff,doBeam,Data::Nt,Data::DoSim);
      }
 #endif
 #ifdef HAVE_CUDA
      if (GPUpredict) {
+      if (beam.elType==ELEM_ALO) {
+        fprintf(stderr,"GPU predict is not supported for this telescope, try CPU only predict\n");
+        exit(1);
+      }
       predict_visibilities_multifreq_withbeam_gpu(iodata.u,iodata.v,iodata.w,iodata.xo,iodata.N,iodata.Nbase,iodata.tilesz,barr,carr,M,iodata.freqs,iodata.Nchan,iodata.deltaf,iodata.deltat,iodata.dec0,
   beam.bfType,beam.b_ra0,beam.b_dec0,beam.p_ra0,beam.p_dec0,iodata.freq0,beam.sx,beam.sy,beam.time_utc,beam.Nelem,beam.xx,beam.yy,beam.zz,&ecoeff,doBeam,Data::Nt,Data::DoSim);
      } else {
@@ -544,11 +548,15 @@ run_fullbatch_calibration(void) {
        predict_visibilities_multifreq(iodata.u,iodata.v,iodata.w,iodata.xo,iodata.N,iodata.Nbase,iodata.tilesz,barr,carr,M,iodata.freqs,iodata.Nchan,iodata.deltaf,iodata.deltat,iodata.dec0,Data::Nt,Data::DoSim);
       } else {
        predict_visibilities_multifreq_withbeam(iodata.u,iodata.v,iodata.w,iodata.xo,iodata.N,iodata.Nbase,iodata.tilesz,barr,carr,M,iodata.freqs,iodata.Nchan,iodata.deltaf,iodata.deltat,iodata.dec0,
-       beam.bfType,beam.b_ra0,beam.b_dec0,beam.p_ra0,beam.p_dec0,iodata.freq0,beam.sx,beam.sy,beam.time_utc,beam.Nelem,beam.xx,beam.yy,beam.zz,&ecoeff,doBeam,Data::Nt,Data::DoSim);
+      beam.bfType,beam.b_ra0,beam.b_dec0,beam.p_ra0,beam.p_dec0,iodata.freq0,beam.sx,beam.sy,beam.time_utc,beam.Nelem,beam.xx,beam.yy,beam.zz,&ecoeff,doBeam,Data::Nt,Data::DoSim);
       }
      }
 #endif
     } else {
+     if (beam.elType==ELEM_ALO) {
+        fprintf(stderr,"GPU predict is not supported for this telescope, try CPU only predict\n");
+        exit(1);
+     }
      /* if solution file is given, read in the solutions and predict */
      read_solutions(sfp,p,carr,iodata.N,M);
 
@@ -734,7 +742,8 @@ run_fullbatch_calibration(void) {
   }
 
   if (doBeam==DOBEAM_FULL||doBeam==DOBEAM_ELEMENT
-      ||doBeam==DOBEAM_FULL_WB||doBeam==DOBEAM_ELEMENT_WB) {
+      ||doBeam==DOBEAM_FULL_WB||doBeam==DOBEAM_ELEMENT_WB
+      ||doBeam==DOBEAM_ALO||doBeam==DOBEAM_ALO_WB) {
    free_elementcoeffs(ecoeff);
   }
   /**********************************************************/
diff --git a/src/MS/main.cpp b/src/MS/main.cpp
index c027514..a25e54e 100644
--- a/src/MS/main.cpp
+++ b/src/MS/main.cpp
@@ -66,6 +66,9 @@ print_help(void) {
    cout << "-z ignore_clusters: if only doing a simulation (with an input solution file using -p option), ignore the cluster ids listed in this file" << endl;
    cout << "-b 0,1 : if 1, solve for each channel: default " <<Data::doChan<< endl;
    cout << "-B 0,1,2,3,4,5,6 : predict "<<DOBEAM_ARRAY<<": array beam, "<<DOBEAM_FULL<<": array+element beam, "<<DOBEAM_ELEMENT<<": element beam, "<<DOBEAM_ARRAY_WB<<": array beam (per channel), "<<DOBEAM_FULL_WB<<": array+element (per channel) "<<DOBEAM_ELEMENT_WB<<": element beam (per channel) : default " <<Data::doBeam<< endl;
+#ifdef HAVE_CSPICE
+   cout << "  For builds with CSPICE (lunar): "<<DOBEAM_ALO<<": ALO element beam, "<<DOBEAM_ALO_WB<<": ALO wideband element beam"<<endl;
+#endif
 #ifdef HAVE_CUDA
    cout << "-E 0,1 : if 1, use GPU for model computing: default " <<Data::GPUpredict<< endl;
 #endif
@@ -147,7 +150,11 @@ ParseCmdLine(int ac, char **av) {
                 break;
             case 'B':
                 doBeam= atoi(optarg);
+#ifdef HAVE_CSPICE
+                if (doBeam>8) { doBeam=DOBEAM_ARRAY; }
+#else
                 if (doBeam>6) { doBeam=DOBEAM_ARRAY; }
+#endif
                 break;
             case 'E':
                 GPUpredict=atoi(optarg);
diff --git a/src/lib/Dirac/Dirac_common.h b/src/lib/Dirac/Dirac_common.h
index 57dad6f..2cc20a5 100644
--- a/src/lib/Dirac/Dirac_common.h
+++ b/src/lib/Dirac/Dirac_common.h
@@ -141,6 +141,14 @@ typedef struct exinfo_shapelet_ {
 #ifndef DOBEAM_ELEMENT_WB
 #define DOBEAM_ELEMENT_WB 6
 #endif
+/* following flags is only to be used for ALO (lunar) simulation,
+ * will need the use of CSPICE */
+#ifndef DOBEAM_ALO
+#define DOBEAM_ALO 7
+#endif
+#ifndef DOBEAM_ALO_WB
+#define DOBEAM_ALO_WB 8
+#endif
 
 typedef struct elementcoff_ {
   int M; /* model order 1,2,3.. */
@@ -261,6 +269,10 @@ typedef struct thread_data_arrayfac_ {
   double *elementgain; /* output : element beam */
   elementcoeff *ecoeff; /* element beam coefficients */
   int dobeam; /* which part of beam to calculate */
+
+#ifdef HAVE_CSPICE
+  pthread_mutex_t *cspice_mutex; /* mutex for CSPICE routines */
+#endif
 } thread_data_arrayfac_t;
 
 
diff --git a/src/lib/Radio/Dirac_radio.h b/src/lib/Radio/Dirac_radio.h
index 0c6db0d..be22ffe 100644
--- a/src/lib/Radio/Dirac_radio.h
+++ b/src/lib/Radio/Dirac_radio.h
@@ -380,6 +380,15 @@ sharmonic_modes(int n0,double *th, double *ph, int Nt, complex double *output);
 #ifdef HAVE_CSPICE
 extern void
 cspice_load_kernels(void);
+
+extern int
+cspice_xyz_to_latlon(double x,double y, double z,double *lon, double *lat, double *alt);
+
+/* calculate element beam, same as element_beam(..) but transforms are in lunar frame
+ * also a mutex is needed as CSPICE is not thread safe
+ */
+extern int
+cspice_element_beam_lunar(double ra, double dec, double f, double f0, int N, double *longitude, double *latitude, double time_jd, elementcoeff *ecoeff, double *elementgain, int wideband, int findex, pthread_mutex_t *mutex);
 #endif /* HAVE_CSPICE */
 
 /****************************** shapelet.c ****************************/
diff --git a/src/lib/Radio/cspice_utils.c b/src/lib/Radio/cspice_utils.c
index 119016d..469df88 100644
--- a/src/lib/Radio/cspice_utils.c
+++ b/src/lib/Radio/cspice_utils.c
@@ -26,6 +26,8 @@
 #include <SpiceUsr.h>
 #include "Dirac_radio.h"
 
+//#define DEBUG
+
 void
 cspice_load_kernels(void) {
   /* get env CSPICE_KERNEL_PATH */
@@ -101,7 +103,142 @@ cspice_load_kernels(void) {
     "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls\n"
     "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/fk/satellites/moon_de440_220930.tf\n"
     "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/moon_pa_de440_200625.bpc\n");
-    fprintf(stderr,"And rerun the program after setting the directory where these kernels stored as CSPICE_KERNEL_PATH\n");
+    fprintf(stderr,"And rerun the program after setting the directory where these kernels stored as CSPICE_KERNEL_PATH in your environment.\n");
     exit(1);
   }
 }
+
+/* x,y,z: selenocentric rectangular coordinates (m)
+ * lon,lat: longitude,latitude (rad)
+ * alt: altitude (m) */
+int
+cspice_xyz_to_latlon(double x,double y, double z,double *lon, double *lat, double *alt) {
+  SpiceDouble radii[3];
+  SpiceInt dim;
+  bodvrd_c("MOON","RADII",3,&dim,radii);
+#ifdef DEBUG
+  printf("MOON RADII %lf %lf %lf\n",radii[0],radii[1],radii[2]);
+#endif
+  /* flattening coefficient */
+  SpiceDouble re=radii[0];
+  SpiceDouble rp=radii[2];
+  SpiceDouble fl=(re-rp)/re;
+  /* rectangular (km) X,Y,Z to planetographic */
+  SpiceDouble rect[3]={x*0.001,y*0.001,z*0.001};
+  recpgr_c("MOON",rect,re,fl,lon,lat,alt);
+
+  (*alt ) *=1000.0; /* km -> m */
+
+  return 0;
+}
+
+
+int
+cspice_element_beam_lunar(double ra, double dec, double f, double f0, int N, double *longitude, double *latitude, double time_jd, elementcoeff *ecoeff, double *elementgain, int wideband, int findex, pthread_mutex_t *mutex) {
+
+  pthread_mutex_lock(mutex);
+  /* ephemeris time from time_jd (days) to (s) */
+  double ep_t0=unitim_c(time_jd,"JED","ET");
+#ifdef DEBUG
+  printf("UTC %le (d) ET %le (s)\n",time_jd,ep_t0);
+#endif
+  /* rectangular coords */
+  SpiceDouble srcrect[3],mtrans[3][3],v2000[3];
+  /* ra,dec to rectangular */
+  radrec_c(1.0,ra*rpd_c(),dec*rpd_c(),v2000);
+  /* precess ep_t0 on lunar frame ME: mean Earth/polar axis, PA: principle axis */
+  pxform_c("J2000","MOON_PA",ep_t0,mtrans);
+  /* rotate v2000 onto lunar frame */
+  mxv_c(mtrans,v2000,srcrect);
+  /* rectangular to lat,long */
+  SpiceDouble s_radius,s_lon,s_lat;
+  reclat_c(srcrect, &s_radius, &s_lon, &s_lat);
+  pthread_mutex_unlock(mutex);
+
+  int ci;
+  double az,el;
+  double theta;
+  /* iterate over stations */
+  for (ci=0; ci<N; ci++) {
+
+  /* find difference between station=lon1,lat1, source=lon2,lat2 */
+  /* Haversine formulae */
+  double d_lon=s_lon-longitude[ci];
+  double d_lat=s_lat-latitude[ci];
+#ifdef DEBUG
+  printf("d_lon %lf %lf %lf\n",s_lon,longitude[ci],d_lon);
+  printf("d_lat %lf %lf %lf\n",s_lat,latitude[ci],d_lat);
+#endif
+  double s_dlat_2=sin(d_lat*0.5);
+  double s_dlon_2=sin(d_lon*0.5);
+  double a=s_dlat_2*s_dlat_2+cos(latitude[ci])*cos(s_lat)*s_dlon_2*s_dlon_2;
+#ifdef DEBUG
+  printf("a=%lf\n",a);
+#endif
+  double a_2=(a>0.0?sqrt(a):1.0);
+  /* great circle distance in rad */
+  double c=2.0*asin((a_2>1.0?1.0:a_2));
+
+  /* azimuth angle, not valid for poles */
+  az=acos( (sin(s_lat)-sin(latitude[ci])*cos(c))/(sin(c)*cos(latitude[ci])) );
+  if (sin(d_lon)>0.0) { az= 2.0*M_PI-az; }
+
+  /* limit latitude difference to [-pi/2,pi/2] */
+  if (fabs(d_lat)<=M_PI_2) {
+  /* elevation = pi/2 - (latitude difference) */
+  el=M_PI_2-fabs(d_lat);
+  } else {
+    el=-1.0; /* invalid */
+  }
+#ifdef DEBUG
+  printf("azimuth %lf elevation %lf\n",az,el);
+#endif
+
+   /* transform : theta = 90-el, phi=-az? 45 only needed for element beam */
+   theta=M_PI_2-el;
+
+   if (el>=0.0) {
+   /* element beam EJones */
+   /* evaluate on r=(zenith angle) 0..pi/2, theta=azimuth grid 0..2pi */
+   /* real data r<- gamma=pi/2-elevation =theta from above code
+    * theta <- beta=azimuth-pi/4  for XX, -pi/2 for YY
+      E = [ Etheta(gamma,beta) Ephi(gamma,beta);
+            Etheta(gamma,beta+pi/2) Ehpi(gamma,beta+pi/2) ]; */
+   if (!wideband) {
+     elementval evalX=eval_elementcoeffs(theta,az-M_PI_4,ecoeff);
+     elementval evalY=eval_elementcoeffs(theta,az-M_PI_4+M_PI_2,ecoeff);
+     elementgain[8*ci]=creal(evalX.theta);
+     elementgain[8*ci+1]=cimag(evalX.theta);
+     elementgain[8*ci+2]=creal(evalX.phi);
+     elementgain[8*ci+3]=cimag(evalX.phi);
+     elementgain[8*ci+4]=creal(evalY.theta);
+     elementgain[8*ci+5]=cimag(evalY.theta);
+     elementgain[8*ci+6]=creal(evalY.phi);
+     elementgain[8*ci+7]=cimag(evalY.phi);
+   } else {
+     elementval evalX=eval_elementcoeffs_wb(theta,az-M_PI_4,ecoeff,findex);
+     elementval evalY=eval_elementcoeffs_wb(theta,az-M_PI_4+M_PI_2,ecoeff,findex);
+     elementgain[8*ci]=creal(evalX.theta);
+     elementgain[8*ci+1]=cimag(evalX.theta);
+     elementgain[8*ci+2]=creal(evalX.phi);
+     elementgain[8*ci+3]=cimag(evalX.phi);
+     elementgain[8*ci+4]=creal(evalY.theta);
+     elementgain[8*ci+5]=cimag(evalY.theta);
+     elementgain[8*ci+6]=creal(evalY.phi);
+     elementgain[8*ci+7]=cimag(evalY.phi);
+   }
+   } else {
+    elementgain[8*ci]=0.0;
+    elementgain[8*ci+1]=0.0;
+    elementgain[8*ci+2]=0.0;
+    elementgain[8*ci+3]=0.0;
+    elementgain[8*ci+4]=0.0;
+    elementgain[8*ci+5]=0.0;
+    elementgain[8*ci+6]=0.0;
+    elementgain[8*ci+7]=0.0;
+   }
+  }
+ 
+
+  return 0;
+}
diff --git a/src/lib/Radio/predict_withbeam.c b/src/lib/Radio/predict_withbeam.c
index 40471af..5a2b90d 100644
--- a/src/lib/Radio/predict_withbeam.c
+++ b/src/lib/Radio/predict_withbeam.c
@@ -190,7 +190,8 @@ precal_threadfn(void *data) {
      }
 
      if (t->dobeam==DOBEAM_ELEMENT || t->dobeam==DOBEAM_FULL
-         ||t->dobeam==DOBEAM_ELEMENT_WB ||t->dobeam==DOBEAM_FULL_WB) {
+         ||t->dobeam==DOBEAM_ELEMENT_WB ||t->dobeam==DOBEAM_FULL_WB
+         ||t->dobeam==DOBEAM_ALO || t->dobeam==DOBEAM_ALO_WB) {
        /* add up terms together with Ejones multiplication */
        for (cn=0; cn<t->carr[cm].N; cn++) {
          complex double *E1=(complex double*)&t->elementbeam[tslot*(t->N*8*t->carr[cm].N)+cn*t->N*8+sta1*8]; /* 8 values */
@@ -410,7 +411,8 @@ precal_threadfn_multifreq(void *data) {
      }
 
      if (t->dobeam==DOBEAM_ELEMENT || t->dobeam==DOBEAM_FULL
-         ||t->dobeam==DOBEAM_ELEMENT_WB ||t->dobeam==DOBEAM_FULL_WB) {
+         ||t->dobeam==DOBEAM_ELEMENT_WB ||t->dobeam==DOBEAM_FULL_WB
+         ||t->dobeam==DOBEAM_ALO || t->dobeam==DOBEAM_ALO_WB) {
        /* add up terms together with Ejones multiplication */
        for (cn=0; cn<t->carr[cm].N; cn++) {
          complex double *E1=(complex double*)&t->elementbeam[tslot*(t->N*8*t->carr[cm].N)+cn*t->N*8+sta1*8]; /* 8 values */
@@ -526,6 +528,19 @@ precalbeam_threadfn(void *data) {
     }
    }
   }
+#ifdef HAVE_CSPICE
+ } else if (t->dobeam==DOBEAM_ALO || t->dobeam==DOBEAM_ALO_WB) {
+    /* iterate over all timeslots */
+    for (ct=0;ct<t->Ntime;ct++) {
+   /* iterate over frequencies */
+   for (cf=0; cf<t->Nf; cf++) {
+   /* iterate over sources */
+  for (cn=t->soff; cn<t->soff+t->Ns; cn++) {
+     cspice_element_beam_lunar(t->carr[cm].ra[cn], t->carr[cm].dec[cn], t->freqs[cf], t->freq0, t->N, t->longitude, t->latitude, t->time_utc[ct], t->ecoeff, &(t->elementgain[ct*(8*t->N*t->carr[cm].N*t->Nf)+cf*(8*t->N*t->carr[cm].N)+cn*8*t->N]),(t->dobeam==DOBEAM_ALO?0:1),cf,t->cspice_mutex);
+    }
+   }
+  }
+#endif
  } else {
   fprintf(stderr,"%s: %d: Invalid mode for beam calculation\n",__FILE__,__LINE__);
  }
@@ -615,7 +630,8 @@ precalculate_coherencies_withbeam(double *u, double *v, double *w, complex doubl
      }
    } 
    if (doBeam==DOBEAM_ELEMENT || doBeam==DOBEAM_FULL
-       ||doBeam==DOBEAM_ELEMENT_WB || doBeam==DOBEAM_FULL_WB) {
+       ||doBeam==DOBEAM_ELEMENT_WB || doBeam==DOBEAM_FULL_WB
+       ||doBeam==DOBEAM_ALO || doBeam==DOBEAM_ALO_WB) {
     /* element beam is common to all stations,
        but az,el might be different, so
        per direction 8 values (2x2 complex): 8N*tilesz, for a cluster with M sources
@@ -789,7 +805,8 @@ precalculate_coherencies_multifreq_withbeam(double *u, double *v, double *w, com
      }
    }
    if (doBeam==DOBEAM_ELEMENT || doBeam==DOBEAM_FULL
-       ||doBeam==DOBEAM_ELEMENT_WB || doBeam==DOBEAM_FULL_WB) {
+       ||doBeam==DOBEAM_ELEMENT_WB || doBeam==DOBEAM_FULL_WB
+       ||doBeam==DOBEAM_ALO || doBeam==DOBEAM_ALO_WB) {
     /* element beam is common to all stations,
        but az,el might be different, so
        per direction 8 values (2x2 complex): 8N*tilesz*Nchan, for a cluster with M sources
@@ -1061,7 +1078,8 @@ visibilities_threadfn_multifreq(void *data) {
      free(tempfr);
 
      if (t->dobeam==DOBEAM_ELEMENT || t->dobeam==DOBEAM_FULL
-         ||t->dobeam==DOBEAM_ELEMENT_WB || t->dobeam==DOBEAM_FULL_WB) {
+         ||t->dobeam==DOBEAM_ELEMENT_WB || t->dobeam==DOBEAM_FULL_WB
+         ||t->dobeam==DOBEAM_ALO || t->dobeam==DOBEAM_ALO_WB) {
       complex double *CO=0;
       if (posix_memalign((void*)&CO,sizeof(complex double),((size_t)t->carr[cm].N*4*sizeof(complex double)))!=0) {
       fprintf(stderr,"%s: %d: No free memory\n",__FILE__,__LINE__);
@@ -1231,6 +1249,9 @@ predict_visibilities_multifreq_withbeam(double *u,double *v,double *w,double *x,
     exit(1);
   }
 
+#ifdef HAVE_CSPICE
+  pthread_mutex_t cspice_mutex=PTHREAD_MUTEX_INITIALIZER;
+#endif
 
   if (add_to_data==SIMUL_ONLY) {
    /* set output column to zero */
@@ -1288,7 +1309,8 @@ predict_visibilities_multifreq_withbeam(double *u,double *v,double *w,double *x,
      }
    }
    if (doBeam==DOBEAM_ELEMENT || doBeam==DOBEAM_FULL
-       ||doBeam==DOBEAM_ELEMENT_WB || doBeam==DOBEAM_FULL_WB) {
+       ||doBeam==DOBEAM_ELEMENT_WB || doBeam==DOBEAM_FULL_WB
+       ||doBeam==DOBEAM_ALO || doBeam==DOBEAM_ALO_WB) {
     /* element beam is common to all stations,
        but az,el might be different, so
        per direction 8 values (2x2 complex): 8N*tilesz*Nchan, for a cluster with M sources
@@ -1336,6 +1358,10 @@ predict_visibilities_multifreq_withbeam(double *u,double *v,double *w,double *x,
      beamdata[nth1].beamgain=beamgain;
      beamdata[nth1].elementgain=elementgain;
      beamdata[nth1].dobeam=doBeam;
+
+#ifdef HAVE_CSPICE
+     beamdata[nth1].cspice_mutex=&cspice_mutex;
+#endif
      pthread_create(&th_array[nth1],&attr,precalbeam_threadfn,(void*)(&beamdata[nth1]));
 
      ci=ci+Nthb;
@@ -1559,7 +1585,8 @@ predict_visibilities_multifreq_withsol_withbeam(double *u,double *v,double *w,do
      }
    }
    if (doBeam==DOBEAM_ELEMENT || doBeam==DOBEAM_FULL
-       ||doBeam==DOBEAM_ELEMENT_WB || doBeam==DOBEAM_FULL_WB) {
+       ||doBeam==DOBEAM_ELEMENT_WB || doBeam==DOBEAM_FULL_WB
+       ||doBeam==DOBEAM_ALO || doBeam==DOBEAM_ALO_WB) {
     /* element beam is common to all stations,
        but az,el might be different, so
        per direction 8 values (2x2 complex): 8N*tilesz*Nchan, for a cluster with M sources
@@ -1846,7 +1873,8 @@ residual_threadfn_multifreq(void *data) {
      }
 
     if (t->dobeam==DOBEAM_ELEMENT || t->dobeam==DOBEAM_FULL
-        ||t->dobeam==DOBEAM_ELEMENT_WB ||t->dobeam==DOBEAM_FULL_WB) {
+        ||t->dobeam==DOBEAM_ELEMENT_WB ||t->dobeam==DOBEAM_FULL_WB
+        ||t->dobeam==DOBEAM_ALO || t->dobeam==DOBEAM_ALO_WB) {
      /* add up terms together with Ejones multiplication */
      for (cn=0; cn<t->carr[cm].N; cn++) {
        complex double *E1=(complex double*)&t->elementbeam[tslot*(t->N*8*t->carr[cm].N*t->Nchan)+cf*(t->N*8*t->carr[cm].N)+cn*t->N*8+sta1*8]; /* 8 values */
@@ -2083,7 +2111,8 @@ calculate_residuals_multifreq_withbeam(double *u,double *v,double *w,double *p,d
      }
    }
    if (doBeam==DOBEAM_ELEMENT || doBeam==DOBEAM_FULL
-       ||doBeam==DOBEAM_ELEMENT_WB || doBeam==DOBEAM_FULL_WB) {
+       ||doBeam==DOBEAM_ELEMENT_WB || doBeam==DOBEAM_FULL_WB
+       ||doBeam==DOBEAM_ALO || doBeam==DOBEAM_ALO_WB) {
     /* element beam is common to all stations,
        but az,el might be different, so
        per direction 8 values (2x2 complex): 8N*tilesz*Nchan, for a cluster with M sources

From 9d2eec562c45a8bfc93a5c903a494cd5c7e8746a Mon Sep 17 00:00:00 2001
From: Sarod Yatawatta <sarod@dop373.nfra.nl>
Date: Sun, 28 Apr 2024 20:16:32 +0200
Subject: [PATCH 12/15] add kernel info

---
 scripts/CSPICE/README.md | 15 ++++++++++++++-
 1 file changed, 14 insertions(+), 1 deletion(-)

diff --git a/scripts/CSPICE/README.md b/scripts/CSPICE/README.md
index 7774dde..70c7170 100644
--- a/scripts/CSPICE/README.md
+++ b/scripts/CSPICE/README.md
@@ -36,4 +36,17 @@ and you are ready to link with sagecal.
 # Linking with sagecal
 Use cmake flag ```-DCSPICE_PREFIX=/full_path_to/cspice``` to point to ```cspice``` directory created after exctracting CSPICE.
 
-za 27 apr 2024  9:23:45 CEST
+
+# Download kernels
+The following kernels are required.
+
+```
+wget https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/pck00011.tpc
+wget https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls
+wget https://naif.jpl.nasa.gov/pub/naif/generic_kernels/fk/satellites/moon_de440_220930.tf
+wget https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/moon_pa_de440_200625.bpc
+```
+
+Download them and save them in a directory, and set your environment variable ```CSPICE_KERNEL_PATH``` to point to this directory before running sagecal.
+
+zo 28 apr 2024 20:16:15 CEST

From 81f1498f480574108aaeb3964c32d086ac63dc0d Mon Sep 17 00:00:00 2001
From: Sarod Yatawatta <sarod@dop373.nfra.nl>
Date: Tue, 30 Apr 2024 12:08:52 +0200
Subject: [PATCH 13/15] change location to load cpice kernels

---
 src/MS/data.cpp           | 4 +---
 src/MS/fullbatch_mode.cpp | 8 ++++++++
 2 files changed, 9 insertions(+), 3 deletions(-)

diff --git a/src/MS/data.cpp b/src/MS/data.cpp
index 261f590..42d1441 100644
--- a/src/MS/data.cpp
+++ b/src/MS/data.cpp
@@ -236,9 +236,7 @@ Data::readAuxData(const char *fname, Data::IOData *data, Data::LBeam *binfo) {
      binfo->elType=(data->freq0<100e6?ELEM_LBA:ELEM_HBA);
    } else if ( !tel.compare("ALO") ) {
      binfo->elType=ELEM_ALO;
-#ifdef HAVE_CSPICE
-     cspice_load_kernels();
-#else
+#ifndef HAVE_CSPICE
      std::cout<<"Warning: telecope "<<tel<<", but CSPICE is not found"<<std::endl;
 #endif
    } else {
diff --git a/src/MS/fullbatch_mode.cpp b/src/MS/fullbatch_mode.cpp
index 9307fd9..3d2b93f 100644
--- a/src/MS/fullbatch_mode.cpp
+++ b/src/MS/fullbatch_mode.cpp
@@ -72,6 +72,14 @@ run_fullbatch_calibration(void) {
     } else if (doBeam==DOBEAM_FULL_WB||doBeam==DOBEAM_ELEMENT_WB||doBeam==DOBEAM_ALO_WB) {
      set_elementcoeffs_wb(beam.elType, iodata.freqs, iodata.Nchan, &ecoeff);
     }
+    if (doBeam==DOBEAM_ALO || doBeam==DOBEAM_ALO_WB) {
+#ifdef HAVE_CSPICE
+      cspice_load_kernels();
+#else
+      std::cout<<"Error: Lunar beam calculation requested by -B option, but CSPICE is not found"<<std::endl;
+      exit(1);
+#endif
+    }
 
 #ifdef HAVE_OPENBLAS
     openblas_set_num_threads(1);//Data::Nt;

From 0ca9b732bb25f1dc84617e17ace5b64b9e5232b7 Mon Sep 17 00:00:00 2001
From: Sarod Yatawatta <sarod@dop373.nfra.nl>
Date: Tue, 30 Apr 2024 13:00:49 +0200
Subject: [PATCH 14/15] update docu

---
 INSTALL.md | 6 +++++-
 README.md  | 2 ++
 2 files changed, 7 insertions(+), 1 deletion(-)

diff --git a/INSTALL.md b/INSTALL.md
index 6c49d13..7f1f2d6 100644
--- a/INSTALL.md
+++ b/INSTALL.md
@@ -1,4 +1,4 @@
-do 15 jun 2023 14:23:00 CEST
+di 30 apr 2024 13:00:36 CEST
 # SAGECal Installation
 
 ## Cmake Build
@@ -31,3 +31,7 @@ SAGECal can use ***libmvec*** vectorized math operations, both in GPU and CPU ve
 ```
 cmake ..  -DCMAKE_CXX_FLAGS='-g -O3 -Wall -ffast-math -lmvec -lm -mavx2' -DCMAKE_C_FLAGS='-g -O3 -Wall -ffast-math -lmvec -lm -mavx2' 
 ```
+
+
+### Linking with CSPICE
+See [linking with CSPICE](https://github.com/nlesc-dirac/sagecal/blob/cspice/scripts/CSPICE/README.md).
diff --git a/README.md b/README.md
index bec8d06..ed8d651 100644
--- a/README.md
+++ b/README.md
@@ -16,6 +16,8 @@
 - Bandpass calibration and unprecedented RFI mitigation with stochastic LBFGS
 - Stochastic calibration for handling data at highest resolution (with federated averaging and consensus optimization)
 - Spectral and spatial regularization of calibration solutions
+- Large scale diffuse sky models via shapelet decomposition
+- Lunar frame interferometric simulation via the SPICE toolkit
 
 Please read INSTALL.md for installation instructions, but 'cmake' should work in most cases. We give a brief guide to use SAGECal here but there is extensive documentation in the links at the end.
 

From f8ca529c94389d85916eb30b1740453e59dbe170 Mon Sep 17 00:00:00 2001
From: Sarod Yatawatta <sarod@dop373.nfra.nl>
Date: Wed, 1 May 2024 09:15:52 +0200
Subject: [PATCH 15/15] fix signs

---
 src/uvwriter/uvwriter.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/uvwriter/uvwriter.cpp b/src/uvwriter/uvwriter.cpp
index aa53ae9..76bff54 100644
--- a/src/uvwriter/uvwriter.cpp
+++ b/src/uvwriter/uvwriter.cpp
@@ -154,8 +154,8 @@ main(int argc, char **argv) {
         SpiceDouble s_radius,s_lon,s_lat;
         reclat_c( srcrect, &s_radius, &s_lon, &s_lat );
         /* [u,v,w]^T=[sinH cosH 0; -sindel*cosH sindel*sinH cosdel; cosdel*cosH -cosdel*sinH sindel] [x y z]^T */
-        double H=s_lon;
-        double del=M_PI_2-s_lat;
+        double H=-s_lon;
+        double del=s_lat;
         rotmat[0][0]=sin(H); rotmat[0][1]=cos(H); rotmat[0][2]=0.0;
         rotmat[1][0]=-sin(del)*cos(H); rotmat[1][1]=sin(del)*sin(H); rotmat[1][2]=cos(del);
         rotmat[2][0]=cos(del)*cos(H); rotmat[2][1]=-cos(del)*sin(H); rotmat[2][2]=sin(del);