Realslide C

realslide.c

/* ********************************* *
*       realslide.c 
*       written: 12 July 2007   
*       updated: 13 July 2007           
*       Jonathan Whitmore
* ********************************* */
#include <stdio.h>
#include <math.h>

int main()
{
int N, i, j, k, l;
double sigma, gprefac;
double niw[1098],nif[1098],nie[1098];   /* phl Non-Iodine wav. and flux */
double iow[1098],iof[1098],ioe[1098];   /* phl Iodine wav. and flux */
double ilw[5258],ilf[5258];             /* Iodine spectra (marcy) wav. and flux */
double con[1098]={0.},new[5258]={0.};
double chivar[61]={0.};

gprefac=1./sqrt(2.*3.14159265358979323846); /* Gaussian prefactor */

FILE *fpni,*fpio, *fmil;        /* phl no iodine, phl iodine, marcy iodine lines */
FILE *fcon, *fish, *fiod;       /* convolved output, iodineshifted, iodine unshifted */
FILE *mtotal, *chisq;           /*  multiplied total, chi squared */
FILE *fcon0,*fcon1,*fcon2,*fcon3,*fcon4;

This next bit changes the scale of the iodine lines:

fmil = fopen("iodinelines.ascii","r");
N=0;
while(fscanf(fmil,"%lf %lf\n",&ilw[N],&ilf[N]) != EOF) N++;

for (i=0;i<=5257;i++)
        {
        ilf[i]=-(ilf[i]-1.60e4)/10.;
        }
/*******        Read in Data Files     ********/
fpni = fopen("noI2_345_PHL957_B67_cropped.ascii","r");
        N=0;
        while (fscanf(fpni,"%lf %lf %lf\n",&niw[N],&nif[N],&nie[N])!= EOF) N++;
fpio = fopen("I2_345_PHL957_B67_cropped.ascii","r");
        N=0;
        while (fscanf(fpio,"%lf %lf %lf\n",&iow[N],&iof[N],&ioe[N]) != EOF) N++;
/*******        Files Ready to output   *********/
fcon = fopen("realconvolution.ascii","w");
mtotal = fopen("realmtotal.ascii","w");
chisq = fopen("realchisq.ascii","w");
fcon1 = fopen("realslideconv1.ascii","w");
fcon2 = fopen("realslideconv2.ascii","w");
fcon3 = fopen("realslideconv3.ascii","w");
fcon4 = fopen("realslideconv4.ascii","w");
/* *********************************** */

This is the part where we do the work
sigma=0.0341;
for (k=0;k<=40;k++)
{
        for (l=0;l<=5257;l++)
        {
                ilw[l]=ilw[l]+0.040000+0.002*k;
        }
                for (i=0;i<=1097;i++)
                {

                        for (j=0;j<=5257;j++)
                        {
                                new[j] = (gprefac/sigma)*ilf[j]*
                                        exp(-((iow[i]-ilw[j])*
                                        (iow[i]-ilw[j])/(2.*sigma*sigma)));
                                con[i] = con[i]+new[j];
                        }

                        if (k==5)
                                fprintf(fcon1,"%lf\t%lf\n",iow[i],con[i]);
                        if (k==8)
                                fprintf(fcon2,"%lf\t%lf\n",iow[i],con[i]);
                        if (k==10)
                                fprintf(fcon3,"%lf\t%lf\n",iow[i],con[i]);
                        if (k==15)
                                fprintf(fcon4,"%lf\t%lf\n",iow[i],con[i]);

                        chivar[k]+=(((nif[i]-iof[i]*1.6)*500.+7.5e4)-con[i])*
                                (((nif[i]-iof[i]*1.6)*500.+7.5e4)-con[i]);

                        con[i]=0;
                }
                printf("%i\n",k);
                fprintf(mtotal,"%lf\t%lf\n",0.040000+0.002*k,chivar[k]);
        for (l=0;l<=5257;l++)
        {
                ilw[l]=ilw[l]-0.04-0.002*k;
        }
}
/* Close files? */
}
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-Share Alike 2.5 License.