#include #include #include #include #include "comphys.h" #include "comphys.c" void sort(unsigned long n, float arr[]); float *scanx,*factorx; float *t; float fx(float x[]); /* Note: Codes 0 = blank image, 1 = star found, but not on slit, 2 = star found and it is on (but may be decentered) the slit. */ int main(int argc, char *argv[]) { float Backx[382] = {235461,235160,234701,235266,234043,234154,234886,233602,233176,234855, 233675,232935,232735,232930,232889,232849,232969,232821,232438,232695, 232539,232605,232184,232331,232406,232690,231730,232186,231565,231736, 231971,232244,231330,232770,231183,232494,231715,231783,231614,231281, 231530,231747,230893,231272,230294,231107,231758,230633,231038,231749, 231188,230847,231380,230410,230759,231026,231740,230935,230674,230543, 230711,230639,230444,230328,230461,229911,230484,230191,229703,229998, 230086,229912,229958,230210,229946,229823,230065,229563,229671,229605, 230011,230004,229175,229406,229836,229482,229575,229126,229531,229237, 229361,228963,229784,228960,230381,228477,229911,229054,229318,229012, 228763,228737,228667,228822,229061,229122,228375,228532,229040,228659, 228396,228529,228558,228421,228274,228366,229155,229424,228704,227996, 228294,228304,228276,228330,228590,228392,228415,228208,227462,227758, 228319,228361,227939,228026,228216,228268,227842,227671,228374,227717, 228061,227682,227703,227893,227883,227623,227735,227351,227628,228085, 227801,227788,227652,226932,227980,227536,228448,227670,227259,228028, 226832,227946,227345,227146,227625,227147,227726,227137,227255,227463, 227326,227663,230613,227171,226864,226904,227435,227227,227239,227455, 227688,227883,227073,226866,227292,227251,227056,226909,227081,226585, 226939,227142,227110,228345,227361,226810,228572,226526,227229,227032, 226862,227351,226874,227089,227013,227095,226922,226851,227329,226890, 226773,227940,226981,226944,227088,226351,227358,226625,227058,226798, 226714,227061,226563,226261,226393,226408,226705,226962,227482,226253, 226643,226411,226293,226724,226576,226724,226459,226683,226429,226740, 226478,226405,226880,226510,226712,226590,226427,226704,226159,226864, 226497,226841,226586,226108,226277,226276,226870,226198,226784,226292, 226912,226296,226358,226525,226183,225863,226502,226469,226788,226359, 225910,226362,227672,226697,226371,226690,226255,226059,226511,226086, 227535,225927,226681,226271,226427,226164,225964,226272,226340,225807, 226102,225994,226240,227177,226716,225996,226668,226199,225779,225877, 226066,226696,226133,225673,226346,226177,225893,227194,225574,226304, 225587,226596,226711,226299,226122,227826,226007,226355,226372,226030, 225952,227796,226684,226356,226743,226058,225957,225745,226143,225615, 226286,226418,225718,226331,225742,225881,228213,225928,226093,226096, 226763,226352,226206,230222,226511,226029,225794,226657,225864,225903, 226274,226055,226178,225706,225670,227656,226283,226242,226515,225627, 225889,226686,225528,226012,226501,225885,225918,225774,226296,226387, 225786,225830,225824,225276,225596,225785,225720,225899,225884,226357, 226484,225657}; float **array; float *stat,*backx; float median,max,maxm1,sumx,sumy,sumM,cx,cy,smax,smax2,xmax2,xmin,xav; int xmax; float medianx,mediany; float smay,ymax,ymax2; float *p,**xi,fret; float xcoord,ycoord; int iter,k; float A,B,C,D,ypos,q; float threshold = 1500.0; int i,j,l,I,J,err,a,b,c,m,n,fwhm; int bx,by; int minfound = 0; char buffer[5000],buf[5000],*tmp,infile[80]; FILE *in,*out,*outy,*bad; backx = Backx - 1; array = matrix(1,382,1,255); scanx = vector(1,382); factorx = vector(1,382); t = vector(1,382); stat = vector(1,97410); p = vector(1,4); xi = matrix(1,4,1,4); out = fopen("scanx.out","w"); strcpy(infile,argv[1]); sprintf(buf,"imlist %s > temp.dat",infile); err = system(buf); for(i=1;i<=382;i++) t[i] = (float)i; in = fopen("temp.dat","r"); if(in == NULL) { printf("Cannot open temp.dat\n"); exit(1); } fgets(buf,4000,in); fgets(buf,4000,in); for(j=255;j>=1;j--) { fgets(buffer,4000,in); tmp = strtok(buffer," "); for(i=1;i<=382;i++) { tmp = strtok(NULL," "); array[i][j] = atof(tmp); } } fclose(in); /* Correct for badpixels by reading badpix.dat file. Skip if it can't find the file */ bad = fopen("badpix.dat","r"); if(bad != NULL) { while(fscanf(bad,"%d,%d",&bx,&by) != EOF) { array[bx][by] = (array[bx-1][by-1] + array[bx][by-1] + array[bx+1][by-1] + array[bx-1][by] + array[bx+1][by] + array[bx-1][by+1] + array[bx][by+1] + array[bx+1][by+1])/8.0; } fclose(bad); } /* Get x-direction scan by adding up columns. Again, this is over the entire chip */ for(i=1;i<=382;i++) scanx[i] = 0.0; for(i=1;i<=382;i++) { scanx[i] = 0.0; for(j=1;j<=255;j++) { scanx[i] += array[i][j]; } } // Background subtraction, x scan for(i=1;i<=382;i++) factorx[i] = scanx[i]/backx[i]; sort(382,factorx); medianx = factorx[191]; for(i=1;i<=382;i++) scanx[i] = scanx[i] - medianx*backx[i]; for(i=1;i<=382;i++) fprintf(out,"%d %f\n",i,scanx[i]); fclose(out); // Find the x-coordinate of the star smax = 0.0; for(j=1;j<=382;j++) { if(scanx[j] > smax) { smax = scanx[j]; xmax = j; } } // No star detected if maximum < threshold if(smax < threshold) { printf("0.00\n"); exit(0); } // If star is in border region, return no star if(xmax < 15 || xmax > 367) { printf("0.00\n"); exit(0); } /* Check -- did we just discover a cosmic ray? Check values +/- 2 pixels from this max. If their average is > 0.5*smax, we have a star. If not, we have a cosmic ray. Fix the cosmic ray and then try again. */ xav = (scanx[xmax-2] + scanx[xmax+2])/2.0; while(xav < 0.5*smax) { scanx[xmax] = 0; smax = 0.0; for(j=1;j<=382;j++) { if(scanx[j] > smax) { smax = scanx[j]; xmax = j; } } xav = (scanx[xmax-2] + scanx[xmax+2])/2.0; } // printf("Xmax = %f smax = %e\n",xmax,smax); /* The Initial Point */ p[1] = A = smax; p[2] = B = xmax; p[3] = C = 10.0; p[4] = D = 0.0; for(i=1;i<=4;i++) { for(j=1;j<=4;j++) { if(i==j) xi[i][j] = 1.0; else xi[i][j] = 0.0; } } /* Invoking Powell */ powell(p,xi,4,0.001,&iter,&fret,fx); xcoord = p[2]; q = p[3]; if(xcoord < 15.0 || xcoord > 367.0) { printf("0.00\n"); exit(0); } printf("%.1f\n",q*0.6); // printf("xcoord = %f\n",p[2]); free_vector(scanx,1,382); free_vector(t,1,382); free_vector(stat,1,97410); free_vector(p,1,4); free_matrix(xi,1,4,1,4); return(0); } float fx(float x[]) { float y1,chi2; int i; chi2 = 0.0; for(i=1;i<=382;i++) { y1 = x[1]*exp(-(t[i]-x[2])*(t[i]-x[2])/(2.0*x[3]*x[3])) + x[4]; chi2 += (y1-scanx[i])*(y1-scanx[i]); } // printf("\nx2 = %f chi2 = %e",x[2],chi2); return(chi2); }