Convolution

Here are two programs that I've written in c to do the calibrations.

This first program calculates a bunch of things, and outputs a giant file of the numbers that it calculates.

The second program looks at the output of the first one, and finds the minimum chi-square for an arbitrary bin size.

/* ********************************* *
 *    convolution.c
 *    written: 10 January 2008 
 *     updated: 24 April 2008
 *    Jonathan Whitmore
 * ********************************* */
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
 
int main(int argc, char *argv[])
{    
    if( argc == 2 )
        printf("Working on order: %s\n", argv[1]);
    else if( argc > 2 )
        printf("Too many arguments supplied.\n");
    else
        printf("One argument expected.\n");
    int N=0, i=0, j=0, k=0, l=0,sizeofbins,numberofshifts;
    int numkecklines=0,nummarcylines=0,sectionlength=0;
    char tablearr[30],noI2[50],I2[50],conmarcy[50],observed[50],chistring[50];
 
    FILE *fpni,*fpio, *fmil;    /* phl no iodine, phl iodine, marcy iodine lines */
    FILE *fcon, *fish, *fiod;     /* convolved output, iodineshifted, iodine unshifted */
    FILE *chisq;            /* chi squared */
    FILE *fcon0,*fcon1,*fcon2,*fcon3,*fcon4; 
    FILE *gnuplot,*shiftsfile;
 
//    printf("Counting the number of lines in the files...");
    numkecklines=0;
    i=0;
    // In case I want to do things with the noI2 data...
    //    sprintf(noI2,"Data/noI2_62.ascii");
    //    fpni = fopen(noI2,"r");
    sprintf(I2,"Data/I2_%s.ascii",argv[1]);
    fpio = fopen(I2,"r");
    N=0;
 
    while(( i = fgetc(fpio)) != EOF)
        if (i == '\n')
            ++numkecklines;
 
    rewind(fpio);
    nummarcylines=0;
    i=0;
    fmil = fopen("Data/splineiodine.ascii","r");
    while(( i = fgetc(fmil)) != EOF)
        if (i == '\n')
            ++nummarcylines;
//    printf("\b\b\b done.\n");
    printf("Number of Keck Lines:  %i\n",numkecklines);
    printf("Number of Marcy Lines: %i\n",nummarcylines);
    /* ************************* */
    double sigma,step=0.,initialshift=0.,twosigsquare=0.;
    double normkeck=0.,sumcon=0.,normratio=0.;
    //    double niw[numkecklines],nif[numkecklines],nie[numkecklines];    /* phl Non-Iodine wav. and flux */
    double iow[numkecklines],iof[numkecklines],ioe[numkecklines];    /* phl Iodine wav. and flux */
    double kdif[numkecklines],differr[numkecklines];        /* Difference in nif and iof*/
    double *ilw,*ilf;
    ilw = (double *) calloc(nummarcylines,sizeof(double));
    ilf = (double *) calloc(nummarcylines,sizeof(double));
 
    printf("Keck data file read: %s\n",I2);
 
    /*******    Read in Data Files     ********/
    N=0;
    while (fscanf(fpio,"%lf %lf %lf\n",&iow[N],&iof[N],&ioe[N])!= EOF) N++;
 
    // Flip Keck lines 
 
    for (i=0;i<numkecklines;i++){
        iof[i]=-iof[i];
    }
 
    double keckoffset=0.;
    double keckarea=0.;
 
    //Iodine Keck Flux Offset
    for (i=0;i<numkecklines;i++){
        keckoffset+=iof[i];
    }
    //printf("Initial Keck Offset: %lf\n",keckoffset);
 
    for (i=0;i<numkecklines;i++){
        iof[i]=  iof[i] - keckoffset/numkecklines;
        if(iof[i]>0){keckarea+=iof[i];}
        if(iof[i]<0){keckarea+=-iof[i];}
    }
    keckoffset=0.;
    for (i=0;i<numkecklines;i++){
        keckoffset+=iof[i];
    }
    //printf("Keck Offset: %lf\nKeck Area: %lf\n",keckoffset,keckarea);
 
    //Find the largest value of normalized flux in the Keck I2
    double bigkeck=0.;
    for (i=0;i<numkecklines;i++){
        if(iof[i] > bigkeck){bigkeck=iof[i];}
    }
//    printf("bigkeck is %lf\n",bigkeck);
 
    /*   Subtracting the lines...    */
    for (i=0;i<numkecklines;i++){
        kdif[i]= iof[i];
    }
 
    /*   Adding errors in quadrature   */
    for (i=0;i<numkecklines;i++){
        differr[i]= sqrt(ioe[i]*ioe[i]);
    }
 
    /* ******     Read in Marcy Iodine **********/
    fmil = fopen("Data/splineiodine.ascii","r");
    N=0;
    while(fscanf(fmil,"%lf %lf\n",&ilw[N],&ilf[N]) != EOF) N++;
    // cut out part of marcy lines
    //    printf("starting the counting process for the iodine lines...\n");
    j=0;
    for ( i = 0; i < nummarcylines; i++ ){
        if ( ilw[i] > iow[0] - 10. && ilw[i] < iow[numkecklines -1] + 10.){
            sectionlength++;
        }
    }
    double section_ilw[sectionlength],section_ilf[sectionlength],marc[sectionlength],convolved[numkecklines];
    for (i=0;i<sectionlength;i++){
        marc[i]=0;
    }
    for (i=0;i<numkecklines;i++){
        convolved[i]=0;
    }
    printf("Size of the iodine section being used: %i\n",sectionlength);
    for ( i = 0; i < nummarcylines; i++ ){
        if ( ilw[i] > iow[0] - 10. && ilw[i] < iow[numkecklines -1] + 10.){
            section_ilw[j] = ilw[i];
            section_ilf[j] = ilf[i];
            j++;
        }
    }
    free(ilw);
    free(ilf);
 
    double bigmarcy=0.;
    for (i = 0 ; i < sectionlength ; i++ ){
        if(section_ilf[i] > bigmarcy){
            bigmarcy=section_ilf[i];
        }
    }
    //printf("bigmarcy is %lf: \n",bigmarcy);
 
    //printf("Changing Iodine lines...");
    for (i = 0 ; i < sectionlength ; i++ ){
        section_ilf[i]=-section_ilf[i]+bigmarcy;
    }
    //printf("\b\b\b done.\n");
 
    //printf("These had better be in increasing order...\n%lf, ",section_ilw[0]);
    //printf("%lf, ",iow[0]);
    //printf("%lf, ",iow[numkecklines-1]);
    //printf("%lf\n",section_ilw[sectionlength-1]);
    sprintf(conmarcy,"Output/marcy_convolved_%s.ascii",argv[1]);
    sprintf(observed,"Output/keck_%s.ascii",argv[1]);
 
    fcon1 = fopen(conmarcy,"w");
    fcon2 = fopen(observed,"w");
 
    /* *********************************** */
    sigma=0.038;
    twosigsquare=(2.*sigma*sigma);
    initialshift=0.010;
    step=0.0001;
 
    shiftsfile = fopen("Parameters/shifts.ascii","r");
    N=0;
    while(fscanf(shiftsfile,"%i\n",&numberofshifts)!= EOF) N++;
 
    printf("\nInitial Shift: %lf\nNumber of shifts: %i\nSize of steps: %lf\n",initialshift,numberofshifts,step);
    printf("Shifted range tested: %lf to %lf\n",initialshift,step*numberofshifts+initialshift);
    double convolvedoffset=0.,convolvedarea=0.,offset=0.,multiply=0.,convolvedareatest=0.;
 
    double **convolvedmatrix;
    convolvedmatrix = (double **)malloc((numkecklines+1)*sizeof(double *));
    for(i=0;i<=numkecklines;i++){
        convolvedmatrix[i] = (double *)malloc((numberofshifts+1)*sizeof(double));
    }
 
    for (l=0;l<sectionlength;l++){
        section_ilw[l]=section_ilw[l]+initialshift;
    }
    //printf("Iodine Lines Initially Shifted, \n");
    for (k=1;k<=numberofshifts;k++){
        if(k%50==0){printf("Iteration: %i; Total Shift = %lf     \n",k,initialshift+step*k);}
        for (l=0;l<sectionlength;l++){
            section_ilw[l]+= step;
        }
        for (i=0;i<numkecklines;i++){
            for (j=0;j<sectionlength;j++){
                if ( (iow[i]-section_ilw[j]) <= 5.*sigma && (iow[i]-section_ilw[j]) >= -5.*sigma ){
                    marc[j] = section_ilf[j]*
                        exp(-((iow[i]-section_ilw[j])*(iow[i]-section_ilw[j]))/(twosigsquare));
                    convolved[i] += marc[j];                    
                }
            }    
        }
        for (i=0;i<numkecklines;i++){
            convolvedoffset+=convolved[i];
        }
        for (i=0;i<numkecklines;i++){
            convolved[i] = convolved[i]-convolvedoffset/numkecklines;
            if(convolved[i]>0){convolvedarea+=convolved[i];}
            if(convolved[i]<0){convolvedarea+=-convolved[i];}
        }
        multiply=keckarea/convolvedarea;
 
        for (i=0;i<numkecklines;i++){
            convolved[i] = convolved[i]*multiply;
        }
 
//        convolvedoffset=0.;
//        convolvedarea=0.;
//        for (i=0;i<numkecklines;i++){
//            convolvedoffset+=convolved[i];
//        }
//        for (i=0;i<numkecklines;i++){
//            if(convolved[i]>0){convolvedarea+=convolved[i];}
//            if(convolved[i]<0){convolvedarea+=-convolved[i];}
//        }
//        printf("Offsets: %lf == %lf\n",convolvedoffset,keckoffset);
//        printf("Areas: %lf == %lf\n",convolvedarea,keckarea);
 
        for (i=0;i<numkecklines;i++){
            fprintf(fcon1,"%lf\t%lf\n",iow[i],convolved[i]);
            convolvedmatrix[i+1][k] = ((kdif[i]-convolved[i])*(kdif[i]-convolved[i])/(differr[i]*differr[i]));
            //    chi squared is with difference squared over sigma^2 
        }
        fprintf(fcon1,"\n\n");
        convolvedareatest=0.;
        convolvedoffset=0.;
        convolvedarea=0.;
        normkeck=0.;
        sumcon=0.;
        for (i=0;i<numkecklines;i++){
            convolved[i]=0;
        }
    }
    for (i=0;i<numkecklines;i++){
        fprintf(fcon2,"%lf\t%lf\t%lf\n",iow[i],kdif[i],differr[i]);
    }
    fprintf(fcon2,"\n\n");
 
    sprintf(chistring,"Output/calculation_%s.ascii",argv[1]);
    chisq = fopen(chistring,"w");
    fprintf(chisq,"%i ",numberofshifts);
    for (i=1;i<=numberofshifts;i++){
        fprintf(chisq,"%lf ",initialshift+i*step);
    }
    fprintf(chisq,"\n");
    for (i=0;i<numkecklines;i++){
        fprintf(chisq,"%lf ",iow[i]);
        for (j=1;j<=numberofshifts;j++){
            fprintf(chisq,"%lf ",convolvedmatrix[i+1][j]);
        }
        fprintf(chisq,"\n");
    }
    for (i=1;i<=numberofshifts;i++){
        convolvedmatrix[0][i]=initialshift+i*step;
    }
    for (i=0;i<numkecklines;i++){
        convolvedmatrix[i+1][0]=iow[i];
    }
    printf("Done with convolving...time to start on the rest.\n");
 
}

Shifting.c

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <error.h>
 
int main(int argc, char *argv[])
{    
    if( argc == 2 )
        printf("Input parameter is: %s\n", argv[1]);
    else if( argc > 2 )
        printf("Too many arguments supplied.\n");
    else
        printf("One argument expected.\n");
 
    int N=0, i=0, j=0, k=0, l=0, numberofshifts=0,numkecklines=0;
    double sizeofbins,input[10];
//Figure out how to read this in
 
    FILE *calculationfile,*binfile,*shiftsfile;    
    binfile = fopen("Parameters/bin.ascii","r");
    N=0;
    while(fscanf(binfile,"%lf\n",&input[N])!= EOF) N++;
 
    printf("Using bin size: %lf\n",input[0]);
    char calcfilestring[100],sortstring[100];
    sprintf(calcfilestring,"Output/calculation_%s.ascii",argv[1]);
    calculationfile = fopen(calcfilestring,"r");
    i=0;
    while(( i = fgetc(calculationfile)) != EOF)
        if (i == '\n')
            ++numkecklines;
 
    rewind(calculationfile);
    shiftsfile = fopen("Parameters/shifts.ascii","r");
    N=0;
    while(fscanf(shiftsfile,"%i\n",&numberofshifts)!= EOF) N++;
 
//numberofshifts=700;
 
printf("Number of lines: %i\nShifts: %i\n",numkecklines,numberofshifts);
 
// Creates memory space for giant calculation matrix.
 
double **convolvedmatrix;
convolvedmatrix = (double **)malloc((numkecklines)*sizeof(double *));
    for(i=0;i<=numkecklines;i++){
    convolvedmatrix[i] = (double *)malloc((numberofshifts)*sizeof(double));
    }
int count;
    N=0;
    for(i=0;i<=numkecklines-1;i++){
        for(j=0;j<=numberofshifts;j++){
            count=fscanf(calculationfile,"%lf ",&convolvedmatrix[i][j]);
            if(count!=1) perror("fscanf conversion\n");
        }
    }
 
sizeofbins=input[0];//remember that the name of the statistics char isn't dynamic
 
int    first = 2;
int bin=0,begin=0,end=0;
double sum[1000][numberofshifts+1],wavelengthofbin[1000];
 
int initialoffset=0;
 
char statistics[200];
 
    sprintf(statistics,"Statistics/sums_%s.ascii",argv[1]);
    sprintf(sortstring,"Statistics/sorthis_%s.ascii",argv[1]);
    FILE *statisticsfile,*outputfile;
    statisticsfile=fopen(statistics,"w");
    outputfile=fopen(sortstring,"w");
 
begin=6;
end=400;
 
for (l=begin;l<=end;l+=13){
first=2;
printf("%i\n",l);
initialoffset=l;
 
    for(i=2;i<numkecklines+1-initialoffset;i++){ 
        if (convolvedmatrix[i+initialoffset][0] - convolvedmatrix[first+initialoffset][0] > sizeofbins){
            bin++;
            for(k=1;k<=numberofshifts+1;k++){
                for(j=first;j<=i;j++){
                    sum[bin][k]+=convolvedmatrix[j+initialoffset][k];
                }
            }
            wavelengthofbin[bin]=(convolvedmatrix[first+initialoffset][0]+convolvedmatrix[i+initialoffset][0])/2.;
            first=i;
 
        }
    }
 
//Here is what the "dynamic" comment is referring to 
    for(i=1;i<=bin;i++){
        fprintf(statisticsfile,"# %lf\n",wavelengthofbin[i]);
        for(j=1;j<=numberofshifts;j++){
            fprintf(statisticsfile,"%lf %lf\n",convolvedmatrix[0][j], sum[i][j]);
        }
        fprintf(statisticsfile,"\n\n");
    }
 
}
 
printf("\n");
 
double min;
int minlocation[bin+1];
 
    for(i=1;i<=bin;i++){
        min = 1000000.;
        for(j=1;j<=numberofshifts;j++){
            if(min > sum[i][j]){
                min = sum[i][j];
                minlocation[i] = j;
            }
        }
 
    }
 
    for(i=1;i<=bin;i++){
        fprintf(outputfile,"%lf  %lf\n",wavelengthofbin[i],convolvedmatrix[0][minlocation[i]]);
    }
 
}
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-Share Alike 2.5 License.