source: Sophya/trunk/SigPredictor/cartelobe.cc@ 1148

// Dominique YVON, CEA/DAPNIA/SPP 02/2000

#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include <string.h>
#include <math.h>
#include <iostream>
//#include "machine.h"
#include "nbconst.h"
#include "cimage.h"
#include "cartelobe.h"


 static NumRecipes NR;  // pour avoir acces au librairies NR

//  Tout cet Objet est programme dans le philosophie des tableaux NR
//  J'essaye d'etre coherent

// Last Index used Value member
long lastXIndex, lastYIndex;

CarteLobe::CarteLobe(char fileRoot[], double frequence, int_4 NbPixelX, int_4 NbPixelY){
  freq=frequence;
  NbPixX=(long)NbPixelX;
  NbPixY=(long)NbPixelY;
  XMin=XMax=YMin=YMax=1.;
  lastXIndex=lastYIndex=1;
  
  pCarte=NR.matrix(1,NbPixX,1,NbPixY);
      // On Reserve l'espace memoire pour stocker la carte. Librairie NR
  pAxeX= NR.vector(1,NbPixX);
  pAxeY= NR.vector(1,NbPixY);  
      // La suite a-t-elle un sens avec les fonction de la NR?
  if ((pCarte==NULL)||(pAxeX==NULL)||(pAxeY==NULL)) {
    cerr<<"erreur memoire dans constructeur pCarte"<<endl;
    cerr<<"espace demandŽ: "<<NbPixX*NbPixY*4<<" Octets"<<endl;
  }

  FILE* pfile;
  char fileCur[150];

  // On lit le tableau des abscisses
  sprintf(fileCur,"%s",fileRoot);
  strcat(fileCur,"_U.txt");
  pfile=fopen(fileCur,"r");
  if (pfile==NULL){
    cerr<< "Erreur a l'ouverture du fichier de donnees :"<<fileCur<< endl;
    exit(-1);
  }

  double PixCur=0.;
  char dumChar;
  for (int_4 i=1; i<(NbPixX+1); i++) {
    fscanf(pfile,"%le%c",&PixCur,&dumChar);
    pAxeX[i]=(r_4)PixCur;
  }

  XMin=XMax=pAxeX[1];
  for (int_4 i=2; i<(NbPixX+1); i++){
    if (XMin>pAxeX[i]) XMin=pAxeX[i];
    if (XMax<pAxeX[i]) XMax=pAxeX[i];
  }
  fclose(pfile);
 
  // On lit le tableau des ordonnees
  sprintf(fileCur,"%s",fileRoot);
  strcat(fileCur,"_V.txt");

  pfile=fopen(fileCur,"r");
  if (pfile==NULL){
    cerr<< "Erreur a l'ouverture du fichier de donnees :"<< fileCur<< endl;
    exit(-1);
  }

  for (int_4 j=1; j<(NbPixY+1); j++) {
    fscanf(pfile,"%le%c",&PixCur,&dumChar);
    pAxeY[j]=(r_4)PixCur;
    for (int_4 i=0; i<(NbPixX-1); i++) {
      fscanf(pfile,"%le%c",&PixCur,&dumChar);
    }
  }

  YMin=YMax=pAxeY[1];
  for (int_4 i=0; i<NbPixY; i++){
    if (YMin>pAxeY[i]) YMin=pAxeY[i];
    if (YMax<pAxeY[i]) YMax=pAxeY[i];
  }
  fclose (pfile);

  // On lit le contenu de la carte 
  sprintf(fileCur,"%s",fileRoot);
  strcat(fileCur,"_PowdB.txt");
  
  pfile=fopen(fileCur,"r");
  if (pfile==NULL){
    cerr<< "Erreur a l'ouverture du fichier de donnees :"<< fileCur<< endl;
   exit(-1);
  }

  for (int j=1; j<(NbPixY+1); j++){
    for(int i=1; i<(NbPixX+1); i++){
      fscanf(pfile,"%le%c",&PixCur,&dumChar);
      pCarte[i][j]=PixCur;
    }
  }
  resolution=CalcResol();
  // C'est fini!
}

CarteLobe::~CarteLobe() {
  NR.free_vector(pAxeX,1,NbPixX);
  NR.free_vector(pAxeY,1,NbPixY);
  NR.free_matrix(pCarte,1,NbPixX,1,NbPixY);
}

double CarteLobe::Value(double XVal, double YVal) const {
        // Value of lobe at coordinates
  float X=(float) XVal;
  float Y=(float) YVal;
  float dBValue=0.;
  float dBError=0.;

  if ((X<=XMin)||(X>=XMax)||(Y<=YMin)||(Y>=YMax)) return 0.;
 
  else {
    unsigned long XIndex=lastXIndex;
    unsigned long YIndex=lastYIndex;
    
    NR.hunt(pAxeX,NbPixX,X,&XIndex);  //  0<X<NbPixX-1
    NR.hunt(pAxeY,NbPixY,Y,&YIndex);  //  0<Y<NbPixY-1

    lastXIndex=XIndex;
    lastYIndex=YIndex;

    // Extrait le soustableau encadrant XIndex et YIndex
    if(XIndex>(NbPixX-2)) XIndex=NbPixX-2;
    if(YIndex>(NbPixY-2)) YIndex=NbPixY-2;
    float ** ptableau= NR.submatrix(pCarte,XIndex,XIndex+2,YIndex,YIndex+2,1,1);
    
    // On interpole sous le soustableau
    NR.polin2(pAxeX+XIndex-1,pAxeY+YIndex-1,ptableau,3,3,X,Y,&dBValue,&dBError);

    // On passe en echelle lineaire 
    return pow(10.,(double)(dBValue/10.));
  }
}

void CarteLobe::FitsVisu(char FitsFile[],long nx,long ny) {
  Image<float>* img;
  img = new Image<float>(nx,ny,1);
  double xstep=(XMax-XMin)/nx;
  double ystep=(YMax-YMin)/ny;
  double yPos=YMin;
  double xPos=XMin;

  for(int i=0; i<nx;i++){
    xPos=XMin;
    for(int j=0; j<ny; j++){
     (*img)(i,j)=Value(xPos,yPos);
      xPos+=xstep;
    }
    yPos+=ystep;
  }
  
  (*img).Save(FitsFile);    // BUGG XXXXXX 
                                                        // Ne sauvegarde pas dans un format FITS pour le moment!!!
  cerr<<"Ecriture fichier :"<<FitsFile<<endl;
  delete img;
}

void CarteLobe::UVToAng(double U, double V, double& alpha, double& beta) const {
  // U=sin alpha cos beta
  // V=sin alpha sin beta
  double sinalpha=sqrt(U*U+V*V);        // Approx des petits angles pour alpha
                                        // Pas pour beta
  if(U>=0) beta=asin(V/sinalpha);
  else beta=Pi-asin(V/sinalpha);
  alpha=asin(sinalpha);
}

double CarteLobe::halfAngAperture(){
  double alCentre=0.; double beCentre=0.; 
  UVToAng(UCentre(),VCentre(),alCentre,beCentre);
  double alBord=0; double beBord=0;
  UVToAng(XMin, YMin, alBord,beBord);
  double CosAng=0.;
  CosAng+=cos(alCentre)*cos(alBord);
  CosAng+=sin(alCentre)*sin(beCentre)*sin(alBord)*sin(beBord);
  CosAng+=sin(alCentre)*cos(beCentre)*sin(alBord)*cos(beBord);
  return acos(CosAng);
}

double CarteLobe::CalcResol(){
  double halfAng=halfAngAperture();
  double NbBin=sqrt(NbPixX*NbPixX+NbPixY*NbPixY);
  return 2*halfAng/NbBin;               // C'est delirament petit
//  on prend
//      return 1.e-3;  // Debug XXXX
}