DOAS.c

Go to the documentation of this file.

#include <windows.h>
#include <stdio.h>
#include <io.h>
#include <fcntl.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <sys/timeb.h>
#include <string.h>
#include <math.h>
#include "mgui.h"
#include "DAS_Spatram.h"
#include "dcl.h"
#include "dil.h"
#include "DOAS.h"
#include "bil.h"

doas DOAS;

int Pix2Wl(int mod, int pos, int spix, int epix, struct doas *pd)
{
    int k, i, c = pos, t=0;

    for(i = 0;i< 18 ;i++)
    {
        if(pos == pd->wlc[i])
        {
                t = i;
        }
    }
    // dl for wl at pix pixref
    pd->lambda[pd->pixref[t]] = (2 * pd->wlref[t] - pd->k1[t]) / (2 - pd->k2[t]); 

    // tot wl at pix pixref
    pd->lambda[pd->pixref[t]] =
        (pd->lambda[pd->pixref[t]]) + (pd->k1[t] - pd->k2[t] * pd->lambda[pd->pixref[t]]) * pd->dhg[t]; 
    

    //Wavelenth series is build for decreasing pixel number
    for(k = pd->pixref[t] - 1 ; k >= spix; k--) 
        pd->lambda[k] = ((2 + pd->k2[t]) * pd->lambda[k + 1] - 2 * pd->k1[t]) / (2 - pd->k2[t]);

    //Wavelenth series is build for increasing pixel number
    for (k = pd->pixref[t] + 1; k <= epix +1; k++) 
        pd->lambda[k] = ((2 - pd->k2[t]) * pd->lambda[k - 1] + 2 * pd->k1[t]) / (2 + pd->k2[t]);

    return 0;
}


 
 
 /*   CalcSpectrumLimits(int mod, int pos, int spx, int epix, int refpix )
     Pixel --> Wavelength
     mod: dummy
     pos: spectral position
     spx: start pixel for wl series
     epx: end pixel for wl series
     refpix: referring pixel
*/
void CalcSpectrumLimits(int mod, int pos, int spix, int epix, int refpix )
{

/*  Da calcoli di Giorgio
2773        0.675070513     -3.42932E-05
3312        0.608620686     -1.22438E-05
3840        0.513266027     9.38545E-06
4358        0.580417561     -8.00615E-06
4866        1.094769302     -1.17723E-04
5363        0.637981086     -2.18565E-05
5850        0.450432353     9.64196E-06
6237        0.926432179     -6.81376E-05
6795        0.509714198     -3.59181E-06
7253        0.961946259     -6.73045E-05
*/

/* Da calcoli di Danbo in DispersioneSPATRAM_EV.xls

  ATTENZIONE LA procedura di linearizzazione utilizza l'equazione:  dl = K1 - K2 * l,
  mentre la procedura di regressione di EXCEL utilizza (giustamente) dl = K1 + K2 * l,
  quindi i valori di k2 in excel hanno il segno "-" mentre qui devono essere inseriti nel
  vettore kk2 con il segno positivo.

2823        dl = 2.35808044 - 0.000628126 * l   
3337        dl = 0.767112208 - 4.67502E-05 * l
3850        dl = 0.824026577 - 5.67659E-05 * l
4358        dl = 3.011942876 - 0.000572933 * l
4861
5332
5814
6282
6734
7201
7677
8195
8622
9023
9405
9758
10066
10358
1.094769302     -1.17723E-04



*/
    
    int k, i, c = pos;
/*  Da calcoli di Giorgio*/
//  double kk1[] = {0.675070513,0.608620686,0.513266027,0.580417561,1.094769302,6,7,8,9,10,11,12,13,14,15,16,17,18};
//  double kk2[] = {-3.42932E-05,-1.22438E-05,9.38545E-06,-8.00615E-06,-1.17723E-04,6,7,8,9,10,11,12,13,14,15,16,17,18};

/* Da calcoli di Danbo*/
//  double kk1[] = {2.35808044,0.767112208,0.824026577,3.011942876,1.094769302,1,1,1,1,1,1,1,1,1,1,1,1,1};
//  double kk2[] = {0.000628126,4.67502E-05,5.67659E-05,0.000572933,1.17723E-04,1,1,1,1,1,1,1,1,1,1,1,1,1};
    
    double lm[1055];

    double k1,k2;

    for(i = 0;i< 18 ;i++)
    {
        if(pos == wl[i])
        {
//          k1 = kk1[i];
//          k2 = kk2[i];
            k1 = DOAS.k1[i];
            k2 = DOAS.k2[i];
        }
    }
/*  
    DOAS.lambda[512] = (2 * pos - k1) / (2 - k2); // dl for wl at pix 512

    DOAS.lambda[512] = (DOAS.lambda[512]) + (k1 - k2 * DOAS.lambda[512]); // tot wl at pix 512
    

    //Wavelenth series is build for decreasing pixel number
    for(k = 511; k >= 0; k--) 
        DOAS.lambda[k] = ((2 + k2) * DOAS.lambda[k + 1] - 2 * k1) / (2 - k2);

    //Wavelenth series is build for increasing pixel number
    for (k = 513; k <= IDX; k++) 
        DOAS.lambda[k] = ((2 - k2) * DOAS.lambda[k - 1] + 2 * k1) / (2 + k2);
*/

    lm[refpix] = (2 * pos - k1) / (2 - k2); // dl for wl at pix pixref

    lm[refpix] = (lm[refpix]) + (k1 - k2 * lm[refpix]); // tot wl at pix pixref
    

    //Wavelenth series is build for decreasing pixel number
    for(k = refpix - 1 ; k >= spix; k--) 
        lm[k] = ((2 + k2) * lm[k + 1] - 2 * k1) / (2 - k2);

    //Wavelenth series is build for increasing pixel number
    for (k = refpix + 1; k <= epix +1; k++) 
        lm[k] = ((2 - k2) * lm[k - 1] + 2 * k1) / (2 + k2);

    for(i = spix; i<=epix; i++)
        DOAS.lambda[i] = lm[i];

}

void FixedStep(void)
{


}



int BuildLogRatio(void)
{

    int     x, y, BlindPix_1;
    FILE *fd;
    float low = 65535;
    float high= -65535;

    IDY = IDY -1;


    BlindPix_1 = 0;
    if (BIL.SPH.ccdwx - (BIL.SPH.ccdex +1) != 0) // se i pixel del sensore sono piu' di quelli in param.ini
        BlindPix_1  = BIL.SPH.ccdwx - (BIL.SPH.ccdex +1) ;

    
//  FixedStep();



    DOAS.IoSmoothed = AllocFloatMat(IDX, IDY);
    DOAS.IsSmoothed = AllocFloatMat(IDX, IDY);

    DOAS.LogRatioMat = AllocFloatMat(IDX, IDY);

    DOAS.LogRatioSmoothed = AllocFloatMat(IDX, IDY);

    DOAS.VarY = AllocFloatMat(IDX, IDY);

    FilterData(DOAS.IoMat, IDX, IDY, 5, BlindPix_1);
    FilterData(BIL.DPLOT.ImatOrig, IDX, IDY, 5, BlindPix_1);
    
    SmoothData(DOAS.IoMat, DOAS.IoSmoothed, IDX, IDY, DOAS.FFTFilter, BlindPix_1);
    SmoothData(BIL.DPLOT.ImatOrig, DOAS.IsSmoothed, IDX, IDY, DOAS.FFTFilter, BlindPix_1);

    for(x = 2; x < IDX - BlindPix_1; x++)
    {
        
        for(y = 0; y < IDY; y++)
        {           
            DOAS.LogRatioMat[y * IDX + x] = (float) log( DOAS.IoMat[y * IDX + x] / ( BIL.DPLOT.ImatOrig[y * IDX + x] * (DOAS.IoSmoothed[y * IDX + x] / DOAS.IsSmoothed[y * IDX + x])));
        }
    }

    SmoothFloatData(DOAS.LogRatioMat, DOAS.LogRatioSmoothed, IDX, IDY, DOAS.FFTFilter, BlindPix_1);

//  for(x = 2; x < IDX - BlindPix_1 + 1; x++)
    for(x = 2; x < IDX - BlindPix_1; x++)
    {
        
        for(y = 0; y < IDY; y++)
        {           
            DOAS.VarY[y * IDX + x] = DOAS.LogRatioMat[y * IDX + x] - DOAS.LogRatioSmoothed[y * IDX + x];
        }
    }


    fd = fopen("TestLogRatio.dat", "w");
    if(fd < 0) return 1;


//  for(x = 2; x < IDX - BlindPix_1 + 1; x++)
    for(x = 2; x < IDX - BlindPix_1; x++)
    {
        
        fprintf(fd, " %9.4lf", DOAS.lambda[x]);
        for(y = 0; y < IDY; y++)
        {   
                
            
//          DOAS.LogRatioMat[y * IDX + x] = log( DOAS.IoMat[y * IDX + x] / ( BIL.DPLOT.ImatOrig[y * IDX + x] * (DOAS.IoSmoothed[y * IDX + x] / DOAS.IsSmoothed[y * IDX + x])));
            fprintf(fd, " %9.4lf", DOAS.VarY[y * IDX + x]);
        }
        fprintf(fd, "\n", NULL);
    }
    fprintf(fd, "\n", NULL);
    fclose(fd);

    



//  for(x = 2; x < IDX - BlindPix_1 + 1; x++)
    for(x = 2; x < IDX - BlindPix_1; x++)
    {
        for(y = 0; y < IDY; y++)
        {
                if(DOAS.VarY[y * IDX + x] < low) low = DOAS.VarY[y * IDX + x];
                if(DOAS.VarY[y * IDX + x] > high) high = DOAS.VarY[y * IDX + x];
        }
    }


//  for(x = 2; x < IDX - BlindPix_1 + 1; x++)
    for(x = 2; x < IDX - BlindPix_1; x++)
    {
        
        for(y = 0; y < IDY; y++)
        {           
            BIL.DPLOT.ImatOrig[y * IDX + x] = (unsigned short)((DOAS.VarY[y * IDX + x] - low) * 100000);
        }
    }
    
    

    

    
//  SwapMatrix((unsigned short *) DOAS.VarY, BIL.DPLOT.ImatOrig, IDX, IDY);
    GrafoIMG(ON_OFFLINEPLOT);

/*
    DOAS.IoSmoothed = DeallocFloatMat(DOAS.IoSmoothed);
    DOAS.IsSmoothed = DeallocFloatMat(DOAS.IsSmoothed);

    DOAS.LogRatioMat = DeallocFloatMat(DOAS.LogRatioMat);

    DOAS.LogRatioSmoothed = DeallocFloatMat(DOAS.LogRatioSmoothed);

    DOAS.VarY = DeallocFloatMat(DOAS.VarY);
*/
    return 0;

}