Commit 5b438c75 authored by Ruben Farinelli's avatar Ruben Farinelli
Browse files

Added the possibility of computing spectropolarimetric results with... 

Added the possibility of computing spectropolarimetric results with estrapolation to lower energy than that defined in the FITS file
parent 679f0bd7

Makefile

+2 −2
Original line number Diff line number Diff line
@@ -22,10 +22,10 @@ endif 
# Sorgenti

SRC = main_program.c hcubature.c restframe_quantities.c polarization_angle.c \

      numpar.c stringpar.c strrev.c interpolate_functions.c trig_funct.c \

      read_file.c sum_values.c read_spectropolarimetry.c array_binning.c

      read_file.c sum_values.c read_spectropolarimetry.c array_binning.c bilinear_extrapolation.c



# Opzioni compilazione

CC = gcc

CC = gcc -g

CFLAGS = -O3

LDFLAGS = -L ${LIB_CFITSIO} -L ${LIB_GSL}  -lcfitsio   -lm   -lgsl   -lgslcblas

HEADERS= -I$(HEADERS_GSL) -I ${HEADERS_CFITSIO}

array_binning.c

+32 −31
Original line number Diff line number Diff line 


#include "headers.h"





void make_array_u(double array_u[], long nx, double start, double end)

{

    if (nx < 2) {

        if (nx == 1) array_u[0] = start;

  if (nx < 2)

  {

    if (nx == 1)

      array_u[0] = start;

    return;

  }

  double step = (end - start) / nx;
@@ -16,10 +17,8 @@ void make_array_u(double array_u[], long nx, double start, double end) 
  }

}





void make_array_energy(double array_energy[], double delta_ene[], double emin, double emax, int nstep)

{

	

  double ymin = log(emin / 511.0);

  double ymax = log(emax / 511.0);

  double dy = (ymax - ymin) / nstep;
@@ -27,13 +26,15 @@ void make_array_energy(double array_energy[], double delta_ene[], double emin, d 
  // Bordi dei bin (nstep + 1 elementi)

  double edges[nstep + 1];



    for (int ii = 0; ii <= nstep; ii++) {

  for (int ii = 0; ii <= nstep; ii++)

  {

    double y_edge = ymin + ii * dy;

    edges[ii] = exp(y_edge) * 511.0;

  }



  // Centri dei bin e larghezza

    for (int ii = 0; ii < nstep; ii++) {

  for (int ii = 0; ii < nstep; ii++)

  {

    array_energy[ii] = 0.5 * (edges[ii] + edges[ii + 1]);

    delta_ene[ii] = edges[ii + 1] - edges[ii];

  }

bilinear_extrapolation.c

0 → 100644
+192 −0
Original line number Diff line number Diff line 
#include "headers.h"



int find_index(double val, double* array, int n)

{

  if (val <= array[0])

    return 0;

  if (val >= array[n - 1])

    return n - 2;



  for (int i = 0; i < n - 1; i++)

  {

    if (array[i] <= val && val <= array[i + 1])

      return i;

  }

  return n - 2; // fallback

}





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



double interp_log_energy_linear_angle(double x, double y, int flag)

{

  double Q11, Q12, Q21, Q22;



  int iu = find_index(y, array_u, nstep_u);

  int ie = find_index(x, array_ene_fits, nstep_ene);



  int all_pos, all_neg;

  

  double LQ1, LQ2, fxy;



 /* printf("Energy %5.4f Angle %5.4f\n", x, y);

  printf("Idx values iu %d ie %d\n", iu, ie);

*/

 

  // coordinate dei punti

  double x1 = array_ene_fits[ie];

  double x2 = array_ene_fits[ie + 1];

  double y1 = array_u[iu];

  double y2 = array_u[iu + 1];



  // valori ai vertici



  if (flag == STOKES_I)

  {

    Q11 = I_energy_angle[iu][ie];         // iu, ie

    Q12 = I_energy_angle[iu][ie + 1];     // iu, ie+1

    Q21 = I_energy_angle[iu + 1][ie];     // iu+1, ie

    Q22 = I_energy_angle[iu + 1][ie + 1]; // iu+1, ie+1

  }

  else

  {

    Q11 = Q_energy_angle[iu][ie];

    Q12 = Q_energy_angle[iu][ie + 1];

    Q21 = Q_energy_angle[iu + 1][ie];

    Q22 = Q_energy_angle[iu + 1][ie + 1];

  }



  all_pos = (Q11 > 0 && Q12 > 0 && Q21 > 0 && Q22 > 0);

  all_neg = (Q11 < 0 && Q12 < 0 && Q21 < 0 && Q22 < 0);

  

  if (all_pos) 

  {

        // log-log positivo

        LQ1 = log(Q11) + (log(Q12)-log(Q11))*(log(x)-log(x1))/(log(x2)-log(x1));

        LQ2 = log(Q21) + (log(Q22)-log(Q21))*(log(x)-log(x1))/(log(x2)-log(x1));

        fxy = LQ1 + (LQ2-LQ1)*(y-y1)/(y2-y1);

        return exp(fxy);

    }

    else if (all_neg) {

        // log-log dei valori assoluti e poi rimettiamo il segno

        LQ1 = log(fabs(Q11)) + (log(fabs(Q12))-log(fabs(Q11)))*(log(x)-log(x1))/(log(x2)-log(x1));

        LQ2 = log(fabs(Q21)) + (log(fabs(Q22))-log(fabs(Q21)))*(log(x)-log(x1))/(log(x2)-log(x1));

        fxy = LQ1 + (LQ2-LQ1)*(y-y1)/(y2-y1);

        return -exp(fxy);

    }

    else {

        // valori misti → interpolazione lineare

        LQ1 = Q11 + (Q12-Q11)*(x-x1)/(x2-x1);

        LQ2 = Q21 + (Q22-Q21)*(x-x1)/(x2-x1);

        fxy = LQ1 + (LQ2-LQ1)*(y-y1)/(y2-y1);

        return fxy;

    }

  

  

  



 

  // printf("Stokes %d CIAOOO %5.4e %5.4e\n", flag, LQ1, LQ2);



  printf("Stokes %d CIAOOO %5.4e %5.4e %5.4e %5.4e\n", flag, Q11, Q12, Q21, Q22);



  // interpolazione lineare lungo l'angolo (y)

  



  // ritorna in scala originale

  return exp(fxy);

}



double bilinear_loglog(double x, double y, int flag)

{

  int iu = find_index(y, array_u, nstep_u);

  int ie = find_index(x, array_ene_fits, nstep_ene);



  double x1 = array_ene_fits[ie];

  double x2 = array_ene_fits[ie + 1];

  double y1 = array_u[iu];

  double y2 = array_u[iu + 1];



  double Q11, Q21, Q12, Q22;



  if (flag == STOKES_I)

  {

    Q11 = I_energy_angle[iu][ie];

    Q21 = I_energy_angle[iu + 1][ie];

    Q12 = I_energy_angle[iu][ie + 1];

    Q22 = I_energy_angle[iu + 1][ie + 1];

  }

  else

  {

    Q11 = Q_energy_angle[iu][ie];

    Q21 = Q_energy_angle[iu + 1][ie];

    Q12 = Q_energy_angle[iu][ie + 1];

    Q22 = Q_energy_angle[iu + 1][ie + 1];

  }



  // logaritmi dei valori

  double LQ11 = log(Q11);

  double LQ21 = log(Q21);

  double LQ12 = log(Q12);

  double LQ22 = log(Q22);



  // interpolazione bilineare in log-log

  double denom = (x2 - x1) * (y2 - y1);

  double Lfxy = LQ11 * (x2 - x) * (y2 - y) / denom + LQ21 * (x2 - x) * (y - y1) / denom +

                LQ12 * (x - x1) * (y2 - y) / denom + LQ22 * (x - x1) * (y - y1) / denom;



  return exp(Lfxy);

}



double bilinear_grid(double x, double y, int flag)

{

  // printf("x=%lf %lf %lf %ld\n",x, array_ene[0], array_ene[nstep_ene-1], nstep_ene );



  double Q11, Q21, Q12, Q22;

  double x1, x2, y1, y2;

  double fxy, denom;

  int iu, ie;



  /* trova gli indici corrispondenti */

  ie = find_index(x, array_ene_fits, nstep_ene); // indice energia (colonna)

  iu = find_index(y, array_u, nstep_u);          // indice u (riga)



  printf("Corresponding idx: iu=%d ie=%d\n", iu, ie);



  /* coordinate dei vertici */

  x1 = array_ene_fits[ie];

  x2 = array_ene_fits[ie + 1];

  y1 = array_u[iu];

  y2 = array_u[iu + 1];



  /* valori ai vertici: ATTENZIONE all'ordine [iu][ie] (rows=u, cols=ene) */

  if (flag == STOKES_I)

  {

    Q11 = I_energy_angle[iu][ie];         // (y1,x1)

    Q21 = I_energy_angle[iu][ie + 1];     // (y1,x2)

    Q12 = I_energy_angle[iu + 1][ie];     // (y2,x1)

    Q22 = I_energy_angle[iu + 1][ie + 1]; // (y2,x2)

  }

  else

  {

    Q11 = Q_energy_angle[iu][ie];

    Q21 = Q_energy_angle[iu][ie + 1];

    Q12 = Q_energy_angle[iu + 1][ie];

    Q22 = Q_energy_angle[iu + 1][ie + 1];

  }



  /* formula bilineare */

  denom = (x2 - x1) * (y2 - y1);



  if (denom == 0.0)

  {

    /* caso degenerato: ritorna il valore al vertice (o fai nearest-neighbour) */

    return Q11;

  }



  fxy = (Q11 * (x2 - x) * (y2 - y) + Q21 * (x - x1) * (y2 - y) + Q12 * (x2 - x) * (y - y1) +

         Q22 * (x - x1) * (y - y1)) /

        denom;



  return fxy;

}

build_model/blcomp_I.fits

−181 MiB (45 MiB)

File changed.

No diff preview for this file type.

View original file
View changed file

build_model/blcomp_Q.fits

−181 MiB (45 MiB)

File changed.

No diff preview for this file type.

View original file
View changed file