File: xsh_eqwidth_lib.c

package info (click to toggle)
cpl-plugin-xshoo 2.8.4%2Bdfsg-3
  • links: PTS, VCS
  • area: main
  • in suites: stretch
  • size: 18,168 kB
  • sloc: ansic: 146,236; sh: 11,433; python: 2,342; makefile: 1,024
file content (893 lines) | stat: -rw-r--r-- 27,416 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

/*-----------------------------------------------------------------------------
                                   Includes
 -----------------------------------------------------------------------------*/

#include "xsh_eqwidth_lib.h"
#include <xsh_msg.h>
#include <gsl/gsl_version.h>

/**@{*/

/*----------------------------------------------------------------------------*/
/**
 * @defgroup xsh_eqwidth_lib
 *
 * TBD
 */
/*----------------------------------------------------------------------------*/


/*----------------------------------------------------------------------------*/
/*
  @brief    get index close to a given wavelenght in a table spectrum
  @param    spec_1d     spectrum table
  @param    lambda      wavelegth
  @return   return index (-1) if not valid????

        
  This function returns a the index of the closest wavelenght greatest than 
  lambda
-------------------------------------------------------------------------------*/
static cpl_size get_index_from_spec(cpl_table* spec_total, double lambda){
    int* error=NULL;
    cpl_size index;
    double cdelta1, crval1;

    cdelta1=cpl_table_get(spec_total, "WAVEL",1, error) - 
                      cpl_table_get(spec_total, "WAVEL",0, error);
    crval1=cpl_table_get(spec_total, "WAVEL",0, error);

    // Obsolete
    //index = (long) (1./cdelta1*lambda-crval1/cdelta1) + 1;
    // Correct formulae for velocity-binned spectra (TO BE USED FOR EFFICIENCY!?)
    //double dv = cdelta1 / crval1;
    //index     = (long)ceil(log(lambda / crval1) / dv);

    // Alternative instructions, using CPL functions
    cpl_table_unselect_all(spec_total);
    cpl_table_or_selected_double(spec_total, "WAVEL", CPL_NOT_GREATER_THAN,
                                 lambda);
    index = (int)cpl_table_count_selected(spec_total);

    if (index <= cpl_table_get_nrow(spec_total))
        {  return index; }
        else { return -1;}
}



static cpl_size get_index_from_spec_alt(cpl_table* spec_total, double lambda){
    int* error=NULL;
    cpl_size index;
    double cdelta1, crval1;

    cdelta1=cpl_table_get(spec_total, "WAVEL",1, error) - 
                      cpl_table_get(spec_total, "WAVEL",0, error);

    crval1=cpl_table_get(spec_total, "WAVEL",0, error);
    index = (long) (1./cdelta1*lambda-crval1/cdelta1) + 1;

// if spectra is binned in velocity 
    double dv = (cdelta1)/crval1;
    double r = 1.+dv;
    index = (long) log(lambda*r/cdelta1)/log(r);

// Python: math.log(ls*(1+dv/c)/ll[0])/math.log(r)

    if (index <= cpl_table_get_nrow(spec_total))
        {  return index; }
        else { return -1;}
}



/*----------------------------------------------------------------------------*/
/**
  @brief    get scale from nm to pixels
  @param    spec_1d     spectrum table
  @param    nm      space scale
  @return   return pixel scale

        
  This function returns a the index of the closest wavelenght greatest than 
  lambda
*/
/*----------------------------------------------------------------------------*/

static cpl_size get_pixel_to_nm_scale(cpl_table* spec_total, double nm){
    int* error=NULL;
    double cdelta1= cpl_table_get(spec_total, "WAVEL",1, error) - 
                      cpl_table_get(spec_total, "WAVEL",0, error);
    return nm/cdelta1;
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Select_local_spec
  @param    spec_1d     spectrum table
  @param    ew_space    parameter for the interval around the line
  @param    line        line center for the spectral region
  @param    spec_1d_out Number of frames in the new frameset
  @return   CPL_ERROR_NONE if OK

        
  This function returns a region of the spectra for the calculations
*/
/*----------------------------------------------------------------------------*/


cpl_error_code select_local_spec(cpl_table* spec_total,
                                  double ew_space,
                                  double lambda,
                                  cpl_table** spec_region)
{
    cpl_errorstate prestate = cpl_errorstate_get();
    int* error = NULL;



    cpl_size index_line = get_index_from_spec(spec_total, lambda);

// 2* space is just a definition...
    cpl_size count = get_pixel_to_nm_scale(spec_total, 2 * ew_space);
    cpl_size begin = index_line - count/2;


    *spec_region = cpl_table_extract(spec_total,
                                     begin,
                                     count);


    if (!cpl_errorstate_is_equal(prestate)) {
        return (int)cpl_error_set_message(cpl_func, cpl_error_get_code(),
               "Unable to Get region of the spectrum");
    }


    return CPL_ERROR_NONE;
}






void find_left_right_continuum_pos(int* x1, int* x2, cpl_table* spec_region,double cont_rejt, double line){



// Interface from CPL tables to C arrays
    int i,n = cpl_table_get_nrow(spec_region);
    double x[n], y[n];
    for(i=0; i < n; i++){
        x[i] = cpl_table_get_double(spec_region,"WAVEL",i,NULL);
        y[i] = cpl_table_get_double(spec_region,"FLUX",i,NULL);
    }


// Code from ARES (could be improved
// It starts from the borders to the inside and check the continuum closest
// location to the line (Confuse but it works)
    int xind1=0,xind2=n-1,hjk;
    double klo=0.01; //Nm instead of Angstroms
    for (hjk=0; hjk < n; hjk++){
        if ((y[hjk] > cont_rejt) && 
            (x[hjk] - (line-klo) > x[xind1] - (line-klo)) && 
            (x[hjk] - (line-klo) < 0))
            xind1=hjk;
        if ((y[hjk] > cont_rejt) && 
            (x[hjk] - (line+klo) < x[xind2] - (line+klo)) && 
            (x[hjk] - (line+klo) > 0) )
            xind2=hjk;
    }
    *x1=xind1;
    *x2=xind2;

}

static cpl_table* esp_spec_deriv(cpl_table* spec_region){

// Interface from CPL tables to C arrays
    int i,n = cpl_table_get_nrow(spec_region);
    double x[n], y[n], dy[n];

    for(i=0; i < n; i++){
        x[i] = cpl_table_get_double(spec_region,"WAVEL",i,NULL);
        y[i] = cpl_table_get_double(spec_region,"FLUX",i,NULL);
    }

    cpl_table* table_deriv = NULL; 	
    table_deriv = cpl_table_duplicate(spec_region);

    deriv(x,y,dy,n);

// Interface from c arrays to CPL Tables:
// Updating the spectral region
    for(i=0; i < n; i++){
        cpl_table_set_double(table_deriv,"FLUX",i,dy[i]);
    }

    return table_deriv;
}


static void esp_spec_smooth(cpl_table* spec_region, int smwidth){

// Interface from CPL tables to C arrays
    int i,n = cpl_table_get_nrow(spec_region);
    double y[n], sy[n];
    for(i=0; i < n; i++){
        y[i] = cpl_table_get_double(spec_region,"FLUX",i,NULL);
    }

    smooth(y, n, smwidth, sy);

// Interface from c arrays to CPL Tables:
// Updating the spectral region
    for(i=0; i < n; i++){
        cpl_table_set_double(spec_region,"FLUX",i,sy[i]);
    }

}




void smooth(double vec[], long n, int w, double svec[]) {
    int i,j;
    double soma;

    if (w%2 != 1)
            w++;
    for (i=0; (i < (w-1)/2);i++)
        svec[i]=vec[i];
    for (i=(w-1)/2;i<n-((w-1)/2);i++) {
        soma=0.;
        for (j=i-((w-1)/2); j<=i+((w-1)/2);j++)
            soma+=vec[j];
        svec[i]=soma/w;
    }
    for (i=n-((w-1)/2); i<n;i++)
        svec[i]=vec[i];
}



static void zeroscenterfind(double iy[], double y[], double dy[], double ddy[], long n, long center[], long *ncenter,double det_line_thres) {
    double maxdy;
    long ctot=0, i, centertot[n];
    int signalc=0, signalc_ant=0;
    if (ddy[0] == abs(ddy[0]))
        signalc=1;
    signalc_ant=signalc;
    maxdy=maxele_vec(dy,n);

    for (i=0; i<n; i++) {
        signalc=0;
        if ( (float) ddy[i] == fabs( (float) ddy[i]) )
            signalc=1;
        // EN: when the signal changes, the local maximum of the 2nd derivative 
        // is below the the noise, and the 3rd derivative is negative enough (due 
        // to noise)
        // EN: for the 0.98 the idea was to have the tree/rejt value, but it does not 
        // work so well. It identifies to many lines in cases of good S/N. This way
        // we only accept identified lines that have a depth of at least 0.98.
        // EN : The 3rd condition was hard codded as iy[i] < 0.98
        if ((signalc != signalc_ant) && 
             (dy[i] > 0.01*maxdy) && 
             (iy[i] < 1. - det_line_thres) && 
             (ddy[i] < -0.1)) {
            centertot[ctot]=i;
            ctot++;
        }
        signalc_ant=signalc;
    }


    if (ctot != 0) {
        *ncenter=ctot;
        for (i=0;i<ctot;i++) 	center[i]=centertot[i];
    } else {
        center[0]=-1;
        *ncenter=0;
    }
}


double maxele_vec(double vec[], long nvec) {
    long i;
    double maxi=vec[1+nvec/20];
    for (i=1+nvec/20; i<nvec-nvec/20; i++)
        maxi = max(maxi,vec[i]);
    return maxi;
}


// NEEDS GSL:::


/*----------------------------------------------------------------------------*/
/**
  @brief    Fit a list of absorption lines with n Gaussian profiles and compute
            the equivalent width of each line from the fit.
  @param    spec_cont_region     spectrum region table (normalized)
  @param    line_table           table with the lines to be fitted
  @param    fit_ngauss_width     Width of the interval used for fitting (nm)???
                                 Currently being used as the guess sigma for 
                                 the gaussian fitting
  @return   CPL_ERROR_NONE if OK

        
  This function detects absorption lines in a spectrum after normalizing it.
*/
/*----------------------------------------------------------------------------*/

cpl_error_code esp_fit_ngauss(cpl_table* spec_cont_region, 
                              cpl_table* line_table,
                              double fit_ngauss_width)
{

    cpl_errorstate prestate = cpl_errorstate_get();


// Interface from CPL tables to C arrays
    int i,ns = cpl_table_get_nrow(spec_cont_region);
    int nl = cpl_table_get_nrow(line_table);

    double xfit[ns], yfit[ns], sigma[ns];

    for(i=0; i < ns; i++){
        xfit[i] = cpl_table_get_double(spec_cont_region,"WAVEL",i,NULL);
        yfit[i] = cpl_table_get_double(spec_cont_region,"FLUX",i,NULL) - 1.0;
//        sigma[i] = cpl_table_get_double(spec_cont_region,"FLUXERR",i,NULL);
//  One good approximation for testing purposes is 1.-cont_rejt
        sigma[i] = 1.- 0.997573;
    }

    int npara=3*nl;
    int igauss=0;
    double acoef[npara], acoef_er[npara];
    for (i=0;i<nl;i++) {
        acoef[3*igauss]=cpl_table_get_double(line_table,"PEAK",i,NULL) - 1.0;
// CAREFUL here with this guess value
        acoef[3*igauss+1]=fit_ngauss_width;
        acoef[3*igauss+2]=cpl_table_get_double(line_table,"WAVEL",i,NULL);
        igauss++;
    }

    int status2;
    
    fitngauss(xfit,yfit,sigma,ns,acoef,acoef_er,npara,&status2);

// We probably need to control here the sucess of the fitting...

// Interface from c arrays to CPL Tables:
// Updating the spectral region
    for(i=0; i < nl; i++){
        cpl_table_set_double(line_table,"PEAK",i,-acoef[3*i]);
        cpl_table_set_double(line_table,"PEAK_ERR",i,acoef_er[3*i]);
        cpl_table_set_double(line_table,"WAVEL",i,acoef[3*i+2]);
        cpl_table_set_double(line_table,"WAVEL_ERR",i,acoef_er[3*i+2]);
// FWHM for the defined Gaussian (ARES): F(X)=Aexp(-Lambda(x-c)^2) => FWHM=2*sqrt(ln(2)/lambda)
        cpl_table_set_double(line_table,"FWHM",i,2.*sqrt(log(2)/acoef[3*i+1]));
        double ind_ew = acoef[3*i]*sqrt(CPL_MATH_PI/acoef[3*i+1])*-1.;
        cpl_table_set_double(line_table,"EW",i,ind_ew*10000.);
        cpl_table_set_double(line_table,"SIGMA",i,sqrt(0.5/(acoef[3*i+1])));
        cpl_table_set_double(line_table,"SIGMA_ERR",i,0.5*CPL_MATH_SQRT1_2*(acoef_er[3*i+1])*pow(acoef[3*i+1],-3.*0.5));        
// using error equation propagation
//    medida_er_square+=medida*medida * ( acoef_er[3*hj]*acoef_er[3*hj]/acoef[3*hj]/acoef[3*hj] + (0.5*0.5*acoef_er[3*hj+1]*acoef_er[3*hj+1]/acoef[3*hj+1]/acoef[3*hj+1]));
//        printf("medida: %15.10lf\n", ind_ew);
        double ind_ew_err_square = ind_ew*ind_ew * ( acoef_er[3*i]*acoef_er[3*i]/acoef[3*i]/acoef[3*i] + (0.5*0.5*acoef_er[3*i+1]*acoef_er[3*i+1]/acoef[3*i+1]/acoef[3*i+1]));
        cpl_table_set_double(line_table,"EW_ERR",i,sqrt(ind_ew_err_square)*10000.);
    }


    if (!cpl_errorstate_is_equal(prestate)) {
        return (int)cpl_error_set_message(cpl_func, cpl_error_get_code(),
               "Unable to Get region of the spectrum");
    }


    return CPL_ERROR_NONE;

}


double check_ew(cpl_table* line_table,
                double line,
                double det_line_resol,
                int* index_line,
                int* n_lines,
                double* ew_error){

    int nl = cpl_table_get_nrow(line_table);

    double ew=0.;
    double ew_er=0.;
	int i;
    
    *index_line = 0;
    *n_lines=0;

	for (i=0; i < nl; i++) {
        if ( fabs( line - cpl_table_get_double(line_table,"WAVEL",i,NULL) ) < det_line_resol ) {
            ew+=cpl_table_get_double(line_table,"EW",i,NULL);
            ew_er+=cpl_table_get_double(line_table,"EW_ERR",i,NULL);
            (*n_lines)++;
            *index_line=i;
        }
    }
    *ew_error=ew_er;

    return ew;
}



static void
poly_fitn(double xvec[], double yvec[], double err[], long n, long ord, double coefs[]) {
    int i, j, k;
    double xi, yi, ei, chisq,xi2;
    gsl_matrix *X, *cov;
    gsl_vector *y, *w, *c;
    ord++;
    X = gsl_matrix_alloc (n, ord);
    y = gsl_vector_alloc (n);
    w = gsl_vector_alloc (n);
    c = gsl_vector_alloc (ord);
    cov = gsl_matrix_alloc (ord, ord);
    for (i = 0; i < n; i++) {
      xi=xvec[i];
      yi=yvec[i];
      ei=err[i];
      for (j = 0; j < ord; j++) {
        xi2=1.0;
        for (k=0; k<j; k++) xi2*=xi;
        gsl_matrix_set (X, i, j, xi2);
        }
      gsl_vector_set (y, i, yi);
      gsl_vector_set (w, i, 1.0/(ei*ei));
    }

    gsl_multifit_linear_workspace * work = gsl_multifit_linear_alloc (n, ord);
    gsl_multifit_wlinear (X, w, y, c, cov, &chisq, work);
    gsl_multifit_linear_free (work);

    #define C(i) (gsl_vector_get(c,(i)))
    #define COV(i,j) (gsl_matrix_get(cov,(i),(j)))

    for (j = 0; j < ord; j++)
        coefs[j]=C(j);
    gsl_vector_free (y);
    gsl_vector_free (w);
    gsl_vector_free (c);
    gsl_matrix_free (X);
    gsl_matrix_free (cov);
    
}

void deriv(double x[], double y[], double dy[], long n) {
    int i;
    gsl_interp_accel *acc = gsl_interp_accel_alloc ();
    gsl_interp *interp = gsl_interp_alloc (gsl_interp_cspline, n);
    //gsl_interp *interp = gsl_interp_alloc (gsl_interp_akima, n);

    gsl_interp_init (interp, x, y, n);
    for (i=0; i<n; i++)
        dy[i]=gsl_interp_eval_deriv (interp, x, y,x[i],acc);
    gsl_interp_free (interp);
    gsl_interp_accel_free (acc);
}


void fitngauss(double t[], double y[], double sigma[], long nvec, 
               double acoef[], double acoef_er[], int para, int *status2)
{
  const gsl_multifit_fdfsolver_type *T;
  gsl_multifit_fdfsolver *s;

  int status;
  size_t i, iter = 0;
  long N=nvec;
  const size_t n = N;
  const size_t p = para;

  gsl_matrix *covar = gsl_matrix_alloc (p, p);
#if defined GSL_MAJOR_VERSION && GSL_MAJOR_VERSION >= 2
  gsl_matrix *J;
#endif
//  double dy[N];
  struct data d = { n, para, t, y, sigma};
  gsl_multifit_function_fdf f;

  double x_init[para];
  for (i=0; i< (size_t) para; i++)
	x_init[i]=acoef[i];

  f.f = &expb_f;
  f.df = &expb_df;
  f.fdf = &expb_fdf;
  f.n = n;
  f.p = p;
  f.params = &d;

  gsl_vector_view x = gsl_vector_view_array (x_init, p);

  T = gsl_multifit_fdfsolver_lmder;

  s = gsl_multifit_fdfsolver_alloc (T, n, p);

  gsl_multifit_fdfsolver_set (s, &f, &x.vector);
  do
    {
      iter++;
      status = gsl_multifit_fdfsolver_iterate (s);

//      int i=0;
      if (status)
        break;

      status = gsl_multifit_test_delta (s->dx, s->x,
                                        1e-6, 1e-6);
    }
  while (status == GSL_CONTINUE && iter < 5000);
#if defined GSL_MAJOR_VERSION && GSL_MAJOR_VERSION >= 2
  J = gsl_matrix_alloc(n, p);
  gsl_multifit_fdfsolver_jac(s, J);
  gsl_multifit_covar (J, 0.0, covar);
  gsl_matrix_free (J);
#else
  gsl_multifit_covar (s->J, 0.0, covar);
#endif

#define FIT(i) gsl_vector_get(s->x, i)
#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))

  {

    double chi = gsl_blas_dnrm2(s->f);
    double dof = n - p;
    double c = GSL_MAX_DBL(1, chi / sqrt(dof));

    for (i=0; i < (size_t) para; i++)
	{
	    acoef[i]=FIT(i);
        acoef_er[i]=c*ERR(i);
	}

  *status2=status;
  *status2=0;   //sometimes we have bad fits but the result is perfectably acceptable
  }
  gsl_multifit_fdfsolver_free (s);
  gsl_matrix_free (covar);
}


// TO BE Located in another file:


/*----------------------------------------------------------------------------*/
/**
  @brief    Create a LINE table with a given row size
  @param    tab  New table
  @param    size Row size
  @return   CPL_ERROR_NONE ff OK

        
  This functions creates a CPL table with structure LINE and a given row size.
 */
/*----------------------------------------------------------------------------*/

cpl_error_code espda_create_line_table(cpl_table **tab,
                                       cpl_size    size) 

{
    cpl_errorstate prestate = cpl_errorstate_get();
    cpl_size       ceil;

/*
 * Table of given row size is created. Columns are those specified for the 
 * structure SPEC.
 */

    *tab = cpl_table_new(size);
    cpl_table_new_column(*tab, "WAVEL",      CPL_TYPE_DOUBLE);
    cpl_table_new_column(*tab, "WAVEL_ERR",  CPL_TYPE_DOUBLE);
    cpl_table_new_column(*tab, "PEAK",       CPL_TYPE_DOUBLE);
    cpl_table_new_column(*tab, "PEAK_ERR",   CPL_TYPE_DOUBLE);
    cpl_table_new_column(*tab, "MU",         CPL_TYPE_DOUBLE);
    cpl_table_new_column(*tab, "MU_ERR",     CPL_TYPE_DOUBLE);
    cpl_table_new_column(*tab, "SIGMA",      CPL_TYPE_DOUBLE);
    cpl_table_new_column(*tab, "SIGMA_ERR",  CPL_TYPE_DOUBLE);
    cpl_table_new_column(*tab, "EW",         CPL_TYPE_DOUBLE);
    cpl_table_new_column(*tab, "EW_ERR",     CPL_TYPE_DOUBLE);
    cpl_table_new_column(*tab, "FWHM",       CPL_TYPE_DOUBLE);
    cpl_table_new_column(*tab, "FWHM_ERR",   CPL_TYPE_DOUBLE);

    if (!cpl_errorstate_is_equal(prestate)) {
        return (int)cpl_error_set_message(cpl_func, cpl_error_get_code(),
               "Unable to create table.");
    }

/*
 * Columns are initialized with default values.
 */

    if (size > 0) {
        ceil = size;
    }
    else {
        ceil = 0;
    }
    cpl_table_fill_column_window_double(*tab, "WAVEL",     0, ceil, -9999.0);
    cpl_table_fill_column_window_double(*tab, "WAVEL_ERR", 0, ceil, -9999.0);
    cpl_table_fill_column_window_double(*tab, "PEAK",      0, ceil, -9999.0);
    cpl_table_fill_column_window_double(*tab, "PEAK_ERR",  0, ceil, -9999.0);
    cpl_table_fill_column_window_double(*tab, "MU",        0, ceil, -9999.0);
    cpl_table_fill_column_window_double(*tab, "MU_ERR",    0, ceil, -9999.0);
    cpl_table_fill_column_window_double(*tab, "SIGMA",     0, ceil, -9999.0);
    cpl_table_fill_column_window_double(*tab, "SIGMA_ERR", 0, ceil, -9999.0);
    cpl_table_fill_column_window_double(*tab, "EW",        0, ceil, -9999.0);
    cpl_table_fill_column_window_double(*tab, "EW_ERR",    0, ceil, -9999.0);
    cpl_table_fill_column_window_double(*tab, "FWHM",      0, ceil, -9999.0);
    cpl_table_fill_column_window_double(*tab, "FWHM_ERR",  0, ceil, -9999.0);

    if (!cpl_errorstate_is_equal(prestate)) {
        return (int)cpl_error_set_message(cpl_func, cpl_error_get_code(),
               "Unable to initialize table.");
    }

/* 
 * The new table is passed along. 
 */

    return CPL_ERROR_NONE;
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Performs a local normalization of the spectrum region
  @param    spec_region     spectrum region table
  @param    cont-rejt   Rejection threshold for the continuum determination.
                        The default value depends on the signal-to­-noise
                        parameter for the interval around the line
  @param    cont-iter   Number of iterations to fit the continuum
  @return   CPL_ERROR_NONE if OK


  This function updates the region of the spectra for the calculations with the
  normalization
*/
/*----------------------------------------------------------------------------*/


cpl_error_code esp_fit_lcont(cpl_table* spec_region,
                             double cont_rejt,
                             int cont_iter){

    cpl_errorstate prestate = cpl_errorstate_get();

    int n = cpl_table_get_nrow(spec_region);


    int order,i,j;
    order=2;
    double coefs[order+1];
    double x[n], y[n], err[n], ynorm[n];
    long nvec;


// Interface from CPL tables to C arrays
    for(i=0; i < n; i++){
        x[i] = cpl_table_get_double(spec_region,"WAVEL",i,NULL);
        y[i] = cpl_table_get_double(spec_region,"FLUX",i,NULL);
//        err[i] = cpl_table_get_double(spec_region,"FLUXERR",i,NULL);
        err[i] = 1.0;
    }


// Calling ARES routines and code (continuum_det5)

    poly_fitn(x,y,err,n,order,coefs);

// Creation of the array with the values of the polynomial
    double xi=1.;
    for(i=0; i < n; i++) {
        ynorm[i]=0.;
        xi=1.;
        for (j=0; j < order + 1; j++) {
            ynorm[i]+=coefs[j]*xi;
            xi*=x[i];
        }
    }

    double vecx[n],vecy[n];
    int jk;
    for (jk=0; jk < cont_iter; jk++) {
        nvec=0;
// Selection of the points for the next fit
        for (i=0; i<n-1; i++) {
//This was to try to avoid cosmic rays, although does not works well always.
// original comment: test were made with 0.01, there should not be any problem to put it larger. I choosed 0.1 in ARES

            if (y[i] > ynorm[i]*cont_rejt && fabs(y[i]-y[i+1]) < 0.1*y[i]){
                vecx[nvec]=x[i];
                vecy[nvec]=y[i];
                nvec++;
            }
        }
        poly_fitn(vecx,vecy,err,nvec,order,coefs);

        for(i=0; i < n; i++) {
            ynorm[i]=0.;
            xi=1.;
            for (j=0;j<order+1;j++) {
                ynorm[i]+=coefs[j]*xi;
                xi*=x[i];
            }
        }
    }

// normalization
    for (i=0; i < n; i++)
        ynorm[i]=y[i]/(coefs[0]+coefs[1]*x[i]+coefs[2]*x[i]*x[i]);

// Interface from c arrays to CPL Tables:
// Updating the spectral region
    for(i=0; i < n; i++){
        cpl_table_set_double(spec_region,"FLUX",i,ynorm[i]);
    }

    if (!cpl_errorstate_is_equal(prestate)) {
        return (int)cpl_error_set_message(cpl_func, cpl_error_get_code(),
               "Unable to Get region of the spectrum");
    }


    return CPL_ERROR_NONE;
}



/*----------------------------------------------------------------------------*/
/**
  @brief    Performs a local normalization of the spectrum region
  @param    spec_cont_region     spectrum region table (normalized)
  @param    det_line_thres       Normalized threshold for line acceptance
                                 The default value is 0.02
  @param    det_line_resol       Minimum accepted separation between two lines
                                 (in Angstrom)
  @param    det_line_smwidth     Number of pixels in the boxcar to be averaged
  @return   CPL_ERROR_NONE if OK


  This function detects absorption lines in a spectrum after normalizing it.
*/
/*----------------------------------------------------------------------------*/

cpl_error_code esp_det_line(cpl_table* spec_cont_region,
                            double det_line_thres,
                            double det_line_resol,
                            int det_line_smwidth,
                            cpl_table** line_table
                            ){

    cpl_errorstate prestate = cpl_errorstate_get();


    cpl_table* spec_cont_region_der1 = NULL;
    cpl_table* spec_cont_region_der2 = NULL;
    cpl_table* spec_cont_region_der3 = NULL;

    spec_cont_region_der1 = esp_spec_deriv(spec_cont_region);

    esp_spec_smooth(spec_cont_region_der1, det_line_smwidth);

    spec_cont_region_der2 = esp_spec_deriv(spec_cont_region_der1);
    esp_spec_smooth(spec_cont_region_der2, det_line_smwidth);

    spec_cont_region_der3 = esp_spec_deriv(spec_cont_region_der2);
    esp_spec_smooth(spec_cont_region_der3, det_line_smwidth);


// Interface from CPL tables to C arrays
    int i,n = cpl_table_get_nrow(spec_cont_region);
    double x[n], iy[n], y[n], dy[n], ddy[n];
    for(i=0; i < n; i++){
        x[i] = cpl_table_get_double(spec_cont_region,"WAVEL",i,NULL);
        iy[i] = cpl_table_get_double(spec_cont_region,"FLUX",i,NULL);
        y[i] = cpl_table_get_double(spec_cont_region_der1,"FLUX",i,NULL);
        dy[i] = cpl_table_get_double(spec_cont_region_der2,"FLUX",i,NULL);
        ddy[i] = cpl_table_get_double(spec_cont_region_der3,"FLUX",i,NULL);
    }

    cpl_table_delete(spec_cont_region_der1);
    cpl_table_delete(spec_cont_region_der2);
    cpl_table_delete(spec_cont_region_der3);

    long ncenter=n, center[n];

//  this obtains the indexes close to the lines center
    zeroscenterfind(iy, y, dy, ddy, n, center, &ncenter,det_line_thres);

//  calculating/interpolating the position of the center of the lines in the spectra:
//  Maybe this is not necessary...

// Only if lines were detected
    if (center[0] != -1 && ncenter != 0) {
        double xlinhas[ncenter], ylinhas[ncenter];
        int i1,i1m;
        double den, diff_y;

        for (i=0; i<ncenter; i++) {
            i1=center[i];
            i1m=i1-1;
            den=1./(x[i1]-x[i1-1]);
            diff_y=(ddy[i1] - ddy[i1m]);
            xlinhas[i]= ( -ddy[i1m] + diff_y*den * x[i1] )  / ( diff_y*den );
            ylinhas[i]= diff_y*den * xlinhas[i] + iy[i1m] - ( iy[i1]- iy[i1m] )*den * x[i1] ;
        }


//RESAMPLING, Negleting "lines" that are to close to each other (noisy lines)
        double xvec2[ncenter], yvec2[ncenter];
        int nvec2,j;
        xvec2[0]=xlinhas[0];
        yvec2[0]=ylinhas[0];
        j=0;
        for(i=1;i<ncenter;i++) {
            if (fabs(xvec2[j]-xlinhas[i]) < det_line_resol ) {
                xvec2[j]=(xvec2[j]+xlinhas[i])*0.5;
                yvec2[j]=(yvec2[j]+ylinhas[i])*0.5;
            } else {
                j++;
                xvec2[j]=xlinhas[i];
                yvec2[j]=ylinhas[i];
            }
        }
        nvec2=j+1;


// Interface C arrays to CPL Table
// Creating and filling the line table

        cpl_ensure_code(espda_create_line_table(line_table,
                                                nvec2)
                        == CPL_ERROR_NONE, cpl_error_get_code());

        for (i=0; i<nvec2; i++) {
            cpl_table_set_double(*line_table,"WAVEL",i,xvec2[i]);
            cpl_table_set_double(*line_table,"PEAK",i,yvec2[i]);
        }
// If no lines are detected then we create a line table without rows
    } else {
        int nvec2=0;
        cpl_ensure_code(espda_create_line_table(line_table,
                                                nvec2)
                        == CPL_ERROR_NONE, cpl_error_get_code());

    }


    if (!cpl_errorstate_is_equal(prestate)) {
        return (int)cpl_error_set_message(cpl_func, cpl_error_get_code(),
               "Unable to Get region of the spectrum");
    }


    return CPL_ERROR_NONE;
}








/**@}*/