#include <stdio.h>
#include <math.h>
#define n 1948
#define nd   7
FILE *fp;
int i1[n],i2[n],i3[n];
float x1[n],x2[n],y;
char a[60];
int d[nd]={105,249,481,725,1041,1621,n};
main(int ac,char **av)
{
  int i,j,k;
  int j0,j1,j2,j3;
  float z1,z2,z3,z4;
  float ava[3],sig[3];
  if(ac<2){
        printf("\n\tcheck ccd?-ccd1 (getfit1 29 24 6) (for 2011)\n\n"); 
        printf("\n\tmyccd1 !  or   myccd1 g\n\n");
        exit(0);
     }
  fp=fopen("/line3/uband-data/cat/fit1.cal1","r");
  k=0;
l10:
  fgets(a,60,fp); if(feof(fp))goto l20;
  if(a[0]=='#')goto l10;
  a[29]=a[34]=a[39]=32;
  sscanf(&a[29],"%d %d %d",&i1[k],&i2[k],&i3[k]);
  x1[k]=x2[k]=0.;
  if(a[49]!=32)sscanf(&a[47],"%f",&x1[k]);
    if(k%100==0)printf("%d %d %d %f\n",i1[k],i2[k],i3[k],x1[k]);
  k++;
  goto l10;
l20:
  fclose(fp);
// ccd1-ccd2, ccd1-ccd3, ccd1-ccd4
  for(i=0;i<k;i++)x2[i]=99.;
  i=0;
l30:
  if(i3[i]!=1){ i++; if(i<k)goto l30; else goto l40; }
  if(x1[i]!=0.){
    if(i3[i+1]==2 && x1[i+1]!=0.)x2[i+1]=x1[i]-x1[i+1];
    if(i3[i+2]==3 && x1[i+2]!=0.)x2[i+2]=x1[i]-x1[i+2];
    if(i3[i+3]==4 && x1[i+3]!=0.)x2[i+3]=x1[i]-x1[i+3];
  }
  i++;
  goto l30;
l40:
  
  j0=0;j1=1;j2=2;j3=3;
  if(av[1][0]!='g')pgbegin_(&j0,"/xw",&j1,&j1,3);
  else { pgbegin_(&j0,"/jg",&j1,&j1,3); pgslw_(&j3);}
  z1=0; z2=k; z3=-.3; z4=0.8;
  pgenv_(&z1,&z2,&z3,&z4,&j0,&j0);
  for(i=0;i<k;i++)if(x2[i]!=99.){
    y=i;
    pgsci_(&i3[i]);
    pgpoint_(&j1,&y,&x2[i],&j1);
  }   
  printf(" 2:\n");
  for(i=j=0;i<k;i++)if(x2[i]!=99. && i3[i]==2)x1[j++]=x2[i];
  sig3(x1,i2,j,&ava[0],&sig[0]);
  printf(" 3:\n");
  for(i=j=0;i<k;i++)if(x2[i]!=99. && i3[i]==3)x1[j++]=x2[i]; 
  sig3(x1,i2,j,&ava[1],&sig[1]);
  printf(" 4:\n");
  for(i=j=0;i<k;i++)if(x2[i]!=99. && i3[i]==4)x1[j++]=x2[i];
  sig3(x1,i2,j,&ava[2],&sig[2]);
  for(i=0;i<nd;i++){
    x1[0]=x1[1]=d[i]-2;
    x2[0]=0.45; x2[1]=0.55;
    pgsci_(&j2); //  if(i==2 || i==8)pgsci_(&j3);
    pgline_(&j2,x1,x2,&j1);
  }  
  z1=300.; pgsci_(&j1);
  z2=0.71; sprintf(a,"Red:   ccd1-ccd2 = %6.3f(mag) rms:%6.3f",ava[0],sig[0]);
  pgtext_(&z1,&z2,a,strlen(a));
  z2=0.66; sprintf(a,"green: ccd1-ccd3 = %6.3f(mag) rms:%6.3f",ava[1],sig[1]);
  pgtext_(&z1,&z2,a,strlen(a));
  z2=0.61; sprintf(a,"blue:  ccd1-ccd4 = %6.3f(mag) rms:%6.3f",ava[2],sig[2]);
  pgtext_(&z1,&z2,a,strlen(a));
  pglabel_("2011: 09_(18,19,20,21,22,23,24)",
   "Mag.","UBAND 4CCD's Flux Ratio",29,4,23);

  pgend_();
}
  
sig3(a,b,m,ava,sig)
float *a,*ava,*sig; int *b,m;          // b is tmp
{
  int i,j,k;
  float x,xx;
  for(i=0;i<m;i++)b[i]=1;
l10:
  x=xx=0.;
  for(i=k=0;i<m;i++)if(b[i]){ x+=a[i]; k++; }
  *ava=x/k; 
  for(i=0;i<m;i++)if(b[i])xx+=(a[i]-*ava)*(a[i]-*ava);
  *sig=sqrt(xx/k);
  printf("m:%4d  av: %6.3f %6.3f\n",k,*ava,*sig);
  x=*sig*3.;
  for(i=j=0;i<m;i++)if(b[i]){
    xx=a[i]-*ava; if(xx<0.)xx=-xx;
    if(xx>x){ b[i]=0; j++; }
  }
  if(j)goto l10;
}