---

APO.h

/*===================================================================
 *
 *    Antenna Positions Optimization
 *
 *   Date: sep 2000
 *   Author: F. Boone
 *   Purpose: Optimize the antenna locations of an interferometer for
 *            a given model distribution of uv-samples.
 *
 *====================================================================*/


#include<stdio.h>
#include<math.h>
#include<time.h>
#include<stdlib.h>
#include<string.h>
#include<iostream.h>                   
#include "cfitsio/fitsio.h"

#define MAXSTRING 50
#define NATOTMAX 500 
#define NDIAMMAX 10 
#define NANTVISIMAX  210000 
#define MAXMAPSIZE   600 // max size of map along one axis
#define MAXDENSITYSIZE 50 // max size of density matrix along one axis
#define NCONFMAX     30 // max # of configurations
#define NPTSMAX      5000 // max # of visibilities for 1 antenna pair
#define NPADSMAX   500// max # of pads
#define PI 3.14159
#define HA2RAD 0.261799 // to transform hour angles into radian
#define RAD2ASEC 206.265e3 // to transform rad into arcsec
#define MINELEVATION 15.0 //minimum source elevation in deg
                          //(should be a free parameter)

/*====================================================================  
 *                            PROTOTYPES
 *====================================================================*/
class Visi; 
class APO;
class Ant;
class Instrume;
class Obs;
class Site;
class Shadow;
class Transf;
class Map;
class Grid;
class Model;
class Configfile;
class Padslist;

extern void checkimplemented(APO *vapo, int napo, 
                             Instrume *pinstr[], Obs *pobs[], Site *psite[], 
                             Grid *pgrid[], Transf *ptransf[],
                             int *sameinstr, int *sameobs, 
                             int *samesite, int *samegrid);
extern int move(APO *vapo, Instrume **pinstr, Obs **pobs, Transf **ptransf,
                Site **psite,
                int napo, int iapo, int iant, float dx, float dy,float dl,
                int sameinstr, int sameobs, int samesite, int samegrid);  
extern void optimize(APO *,int,float,float);
extern void optimize_shad(APO *,int,float,float, float gsh);
extern void optimize(APO *,int,float,float, int conf);
extern void shake(APO *vapo,int napo, float dptacc);
extern void correct(APO *vapo,int napo, float dptacc);
extern void correct(APO *vapo,int napo, float dptacc, int conf);
extern float solveFWHM(float B, float Bmax);
extern float bilinear(float, float, float **);
extern float **matrix(int, int);
extern void delete_matrix(float **);
extern int compi(const void *, const void *);
extern int compf(const void *, const void *);
extern int compvec(const void *v1, const void *v2);
extern int compdist(const void *v1, const void *v2);
extern void printfitserror( int status);

/*====================================================================  
 *                            CLASSES
 *====================================================================*/


class Visi{
  float u,v,r,t,ix,iy,wght;
public:
  inline Visi();
  inline void set(float, float, float);
  inline void set(float, float, float, float, float);
  inline void set(float, float, float, float, float, float, float);
  inline void operator=(Visi);
  inline void operator=(float);
  inline void getijwght(int *, int *, float *);
  inline void getuvwght(float *u, float *v, float *wght);
  inline void sym(Grid *);

  void print(){
    printf("u=%7.1f v=%7.1f r=%7.1f t=%7.4f wght=%5.3f\n",u,v,r,t,wght);
    printf("                   ix=%7.2f iy=%7.2f\n",ix,iy);
  };
};
  
/*====================================================================*/  

class APO{

  long nvisi, nvisitot;
  Grid *pgrid;
  Instrume *pinstr;
  Obs *pobs;
  Site *psite;
  Transf *ptransf;
  float **pweights;
  float **pdensity; 
  float **pgradx;  
  float **pgrady;   
  float **dpt; 
  Visi *pvisi;                 

  inline void getvisi(int);
  inline void initweights(int);
  inline void addallvisi();

public:
  APO();
  ~APO();
  void set(Instrume *, Obs *, Transf *, Site *, Grid *);
  void set(Grid *);
  void reset();
  void correct(float dl);
  void writeuvcov(char* fitsfile, long npix, float du);
  inline void addvisi(int);
  inline void subtractvisi(int);
  inline void getgrad();
  inline int check(int iant, float x, float y,int nconf1,int *iconf1);
  float getstdev();
  void shake(float dptacc);
  inline int move(int, float, float, float);
  void optimize(float, float);
  void optimize(float, float, int);
  void printvisi();
  void printdensity();
  void printforce(int iant, float fa, char *filename);
  void printforce(int iant);
  inline void getforce(float *, float *,int, float);
  inline void getforce_shad(float *, float *,int, float, float gsh);
  Instrume *getpinstr(){return pinstr;};
  Obs *getpobs(){return pobs;};
  Transf *gettransf(){return ptransf;};
  Site *getpsite(){return psite;};
  Grid *getpgrid(){return pgrid;};

};
  
/*====================================================================*/  


class Ant{ 
  float x;
  float y;
  float diam;
  int fixed;
  int moved;
  int nconf;
  int iconf[3];
public:
  Ant(){;};
  void set(float a, float b, int fixed1, float c, int d, int e[]){
    x=a; y=b; fixed=fixed1; diam=c; nconf=d;
    for (int i=0; i<nconf; i++) iconf[i]=e[i]; moved=0;
  };  
  void set(float a, float b){
    x=a; y=b; 
  };  
  void operator=(Ant ant2){
    int i;
    x=ant2.x; y=ant2.y; diam=ant2.diam; nconf=ant2.nconf; 
    for(i=0;i<nconf;i++)iconf[i]=ant2.iconf[i];
  };
  void get(float *p,float *q, int *pfixed, float *r){
    *p=x; *q=y; *pfixed=fixed; *r=diam;
  };
  void get(float *p,float *q){*p=x; *q=y;};
  void get(float *p,float *q, int *pfixed){*p=x; *q=y; *pfixed=fixed;};
  void get(float *p, float *q, int *pnconf1, int *piconf1){
    int i;
    *p=x; *q=y; *pnconf1=nconf;
    for(i=0;i<nconf;i++) *(piconf1+i)=iconf[i];
  };
  void move(float dx, float dy){
    x+=dx; y+=dy; moved=(moved+1) % 2;
  };
  int getmoved(){return moved;};
  int getfixed(){return fixed;};
  float getdiam(){return diam;};
  void resetmoved(){moved=0;};
  void freeant(){fixed=0;};
  inline int check(Site *);
  void fixant(float x1, float y1){x=x1;y=y1;fixed=1;};
  void fixant(){fixed=1;};
  void offset(float xoff, float yoff);
  friend void APO::getvisi(int);  
  void print(){  
      int i;
      printf("x= %7.1f   y=%7.1f   diam=%4.1f   conf:",x,y,diam);
      for(i=0; i<nconf; i++) printf(" %d",iconf[i]);
      printf(" fixed=%d\n",fixed);
  };
};
  
/*====================================================================*/ 


class Instrume{
  Ant *AntList;
  int na, natot, ndiam, nsh, nconf, nmoved, notmoved;
  float dspec; 
  float diamlist[NDIAMMAX];
  float *vdmax;
  float *vdmin;

public:
  Instrume(){
    AntList=NULL; vdmax=NULL; vdmin=NULL, nmoved=0;
  };
  ~Instrume(){
    delete AntList; delete vdmax; delete vdmin;
    AntList=NULL; vdmax=NULL; vdmin=NULL;
  };
  void set(float *, float *, int *fixed1, float *, int, int, int, float, 
           float *, float *);
  void set(float *, float *, int na);
  void setvdmax(float *vdmax1, int nconf1);
  void random(Site *psite, Obs *pobs, Model *pmod);
  void sort();
  void offset(float xoff, float yoff);
  void fixant(int iant){AntList[iant].fixant();};
  inline void get(float *, float *, int *fixed1, float *, int *, 
                  int *, int *, float *, float *, float *);
  void getpositions(float *xx, float *yy);
  void getcenter(float *xc, float *yc);
  inline void getant(float *p, float *q, int *pnconf1, int *piconf1,int iant){
    AntList[iant].get(p, q, pnconf1, piconf1);
  };
  inline Ant getant(int iant){return AntList[iant];};
  void getcorrection(Transf *ptransf, Site *psite, int iant0,
                     int *iant1, float *dx, float *dy);
  void getcorrection(Transf *ptransf, Site *psite, int iant0, int conf,
                     int *iant1, float *dx, float *dy);
  int getnatot(){return natot;};
  int getna(){return na;};
  int getnconf(){return nconf;};
  int getnsh(){return nsh;};
  int getiant(int conf, int i){return conf*(na-nsh)+i;};
  int getmoved(int iant){return AntList[iant].getmoved();};
  int getfixed(int iant){return AntList[iant].getfixed();};
  int getndiam(){return ndiam;};
  float getdiam(int i);
  float getdmax();
  inline float getdptmax(int iant,float x1, float y1, float dx ,float dy,
                        float dl, int nconf1,int iconf1[3], Transf *ptransf);
  int getBSL_FWHM_ObsDur(float obsdurmax, char *filename, int iconf, 
                          float *baseline, float *fwhm, float *obsdur);
  inline int check(int, float, float, int, int *, Transf *);
  inline int check(int, float, float, int, int *);
  inline int check(Transf *);
  inline int check(Site *);
  inline void move(int iant, float dx, float dy);
  inline void resetmoved();
  void freeants();
  //void random(Site *);
  void addants2pads(Padslist *pads, int iconf);
  int loadants2pads(Padslist *pads, int nloadmax, float dist, 
                    float ff,int iconf);
  void print();
  friend APO; 
};
 
  
/*====================================================================*/ 
 

class Obs{
  float dec;
  float dt;
  int nconf;
  float *vhamin;
  float *vintime;
  int *vnpts;
public:
  Obs(){
    vhamin=NULL; vnpts=NULL;
  };
  ~Obs(){
    delete vhamin;
    delete vnpts;
    vhamin=NULL; vnpts=NULL;
  };
  void Obs::operator=(Obs obs2);
  void set(float, float *, float, float *, int);
  void set(float dec1);
  void get(float *, float *, float *, float *, int *);
  void get(float *hamin,float *hamax,int ind){
    *hamin=vhamin[ind]*HA2RAD;
    if(vnpts[ind] > 1) *hamax=(vhamin[ind]+vintime[ind])*HA2RAD;
    else *hamax=*hamin;
  };
  float getdec(){return dec;};
  float getobsdur(int conf){return vintime[conf];}; 
  int getnpts(int conf){return vnpts[conf];};
  void print();
  friend void APO::getvisi(int);  
  friend void APO::initweights(int);
};
/*====================================================================*/  

class Shadow{
  float sinlat, sindec, cosdl,sindl;
  float *pcosha;
  float *psinha;
  int npts;
 public:
  Shadow();
  ~Shadow();
  void set(float dec,float lat, float elevmin);
  inline float getruvmin(float blx, float bly);
};
/*====================================================================*/  

class Transf{
  Obs *pobs;
  Site *psite;
  int nconf;
  float coslat, sinlat, cosdec, sindec, cosdl,sindl, cosd_l;
  float **pcosha;
  float **psinha;
  int *vnpts;
  int shadowing;
  Shadow shad;
 public:
  Transf();
  ~Transf();
  void set(Obs *pobs, Site *psite);
  void setshadow(float dec);
  inline void xy2uv(float blx, float bly, int conf, int isample, 
             float *pu, float *pv);
  inline void uv2xy(float u, float v, int conf, int isample, 
             float *px, float *py);
  inline float getruvmax(float blx, float bly, int conf);
  inline float getruvmin(float blx, float bly, int conf);
  float getcosha(int isample,int conf){return pcosha[isample][conf];};
  void get(float *coslat1, float *sinlat1, float *cosdec1, float *sindec1, 
                 float *cosdl1, float *sindl1, float *cosd_l1);
};
/*====================================================================*/  

class Map{  

  long naxis1, naxis2;
  float xcenter, ycenter;
  float cdelt1, cdelt2;
  float **parray;
  char filename[MAXSTRING];

public:
  Map(){parray=NULL;};
  ~Map(){if(parray!=NULL) delete_matrix(parray);};
  inline int check(float x, float y);
  void set(char *filename1, float xcenter1,float ycenter1,
           float cdelt11, float cdelt21);
  void get(char *filename1,float *xcenter1,float *ycenter1,
           float *cdelt11,float *cdelt21){
    strcpy(filename1,filename);
    *xcenter1=xcenter; *ycenter1=ycenter;
    *cdelt11=cdelt1; *cdelt21=cdelt2;
  };
};
  
/*====================================================================*/  

class Site{

  Map ground;
  int nomap;
  float lat;
  float tau;
  float rmin, rmax;     // min & max distance possible from center of site
  float elevmin;      // minimum possible elevation in deg
public:
  Site(){;};
  void set(char *filename1, float xcenter1,float ycenter1, 
           float cdelt11, float cdelt21,
           float lat1, float tau1,float rmin1,float rmax1);
  void set(float lat1, float tau1, float rmin1,float rmax1){
    lat=lat1; tau=tau1; nomap=1; rmin=rmin1; rmax=rmax1;
    printf("nomap=%d\n",nomap);  
    elevmin=MINELEVATION;
  };
  void get(float *lat1, float *tau1, 
           float *rmin1, float *rmax1, int *nomap1){
    printf("nomap=%d\n",nomap);
    *lat1=lat; *tau1=tau; *nomap1=nomap; *rmin1=rmin; *rmax1=rmax;
  };
  void get(char *filename1,float *xcenter1,float *ycenter1,
           float *cdelt11,float *cdelt21){
    ground.get(filename1, xcenter1, ycenter1, cdelt11, cdelt21);
  };
  inline int check(float, float);
  float getlat(){return lat;};
  float getelevmin(){return elevmin;};
  void print();
  friend void APO::getvisi(int); 
  friend void APO::initweights(int);
};

/*====================================================================*/  

class Model{
  float fwhm, sigsq, dmax, dmin, norm;
public:
  Model(){;};
  void set(float, float, float);
  float integral(float r0,float r1,float daz){
    return (norm*daz*sigsq*(exp(-r0*r0/(2*sigsq))
                         -exp(-r1*r1/(2*sigsq))));
  };
  void setfwhm(float fwhm1){fwhm=fwhm1;};
  float getfwhm(){return fwhm;};
  void get(float *fwhm1, float *dminmod1,float *dmaxmod1){
    *fwhm1=fwhm; *dminmod1=dmin; *dmaxmod1=dmax;
  };
  float getdmax(){return dmax;};
  float getdmin(){return dmin;};
  float getnextr(float r0,float integr, float daz){
    return sqrt(-2*sigsq*log(exp(-r0*r0/(2*sigsq))-integr/(daz*norm*sigsq)));
  };
};

/*====================================================================*/ 

class Grid{
  int size;
  float integr;  // value of the integral of the model over 1 bin
  float daz;
  float radius[MAXDENSITYSIZE+1];
  float az[2*MAXDENSITYSIZE+1];
  Model *mod;
public:
  Grid(){;};
  void set(int, float, float, Model *);
  void setrmax(float rmax);
  inline void getixiy(float *, float *, float, float, float, float);
  inline float getdr(float iy);
  int getsize(){return size;};
  float getdaz(){return daz;};
  float getintegr(){return integr;};
  float getfwhm(){return mod->getfwhm();};
  void print();
  void print(char *name);
};

  
/*====================================================================*/  

class Configfile{
  FILE *fp;
  char name[40];
  void read(int *tmp);
  void read(long *tmp);
  void read(float *tmp);
  void read(char *s);
  void read(float *pv,int n);
  void read(int *pv,int n);
public:
  Configfile(){
    fp=NULL;
  };
  void operator=(char *name1){
    strcpy(name,name1);
  };
  void read(Instrume *, Site *, Obs *, Model *);
  void write(Instrume *, Site *, Obs *, Model *);
};

  
/*====================================================================*/  

class Padslist{
  char filename[MAXSTRING];
  long npads;
  int nconf;
  float **pos;
  int **conf;
  
public:
  Padslist();
  ~Padslist();
  void setnconf(int nconf1){nconf=nconf1;};
  char *getfilename(){return filename;};
  void operator=(char *filename1){strcpy(filename,filename1);};
  float getdist2pad(float x, float y, int iconf, 
                    float *xpad, float *ypad, int *ipad);
  int load(int ipad, int iconf);
 // int load(float x, float y, float dist, int iconf, float *xpad, float *ypad);
//  long readnpads();
  void add(float x, float y, int iconf);
  void read();
  void write();
  void write(char *filename2);
  void write2(char *filename2);
  void offset(float dx, float dy);
  void print();
};

---