/* zusammenskal.c     Programm haengt blockweise gemessene Daten aus HASYLAB
aneinander zu einer DAT Datei fuer SIMREF und benutzt ausserdem einen vom 
Benutzer erfragten und in der Datei mit den gemessenen Winkelbereichen immer 
gleichartig vorkommenden Referenzwinkelbereich, um die Intensitaetsdriften am 
HASYLAB durch Korrekturen an den Intensitaeten auszugleichen.
Es wird die Datei nach diesen Referenzmessungen durchsucht, und jeweils
die Summe der Zaehlungen mit der Summe des ersten Bereichs verglichen,
um den Korrekturfaktor zu finden.

Uebersetzen mit: cc -lm zusammen.c  (sonst findet er die Mathe-Funktionen nicht)

Eingabe: * Dateinamen der Eingabedatei im HASYLAB Datenformat.
         Sie enthaelt Abfolge von Bloecken in der richtigen 
         Reihenfolge zu groesseren 2 theta Winkeln hin.         
         Wenn die Bloecke keine gleiche 2-theta Schrittweite besitzen,
         koennen sie dennoch eingeordnet werden!
         Die gleiche Abfolge von Winkelbereichen kann auch spaeter in
         der Datei wiederholt werden, was zum 'mergen' der Daten benutzt wird
         und die Genauigkeit der Daten steigert.
         Beliebige 2theta-Werte in aufsteigender Reihenfolge koennen 
         vorkommen, eine feste Schrittweite muss nicht eingehalten werden.
Ausgabe: * Eine Datei mit nur einem Block, mit den aneinandergehaengten
         Messwerten (xxxxx.zus) im DAT-Format: 
               2theta, Messwert, sigma(Messw), Untergrundwert.

         * Eine Datei mit den EXCLUDED REGIONS (xxx.log), dass heisst
         den Winkelbereichen, die garnicht gemessen wurden.
  
Erstellt durch K. Burger, Okt. 1994 (zusammen.c), modifiziert fur DAT
statt DUO Dateien Sept. 1996, modifiziert fuer automatische interne
Intensitaetskorrektur bzw. -skalierung Mai 1997.

*/

/* >>>>>>>>>>>>>>   I N C L U D E -- D A T E I E N   <<<<<<<<<<<<<<<<<< */

#include "/usr/include/stdio.h"
#include "/usr/include/stdlib.h"
#include "/usr/include/math.h"
#include "/usr/include/string.h"
#define MAXWERTE 2000
#define MAXGESWERTE 30000
#define MAXBLOCK 200

int redeviel=1;     /* nur zum Testen */
float banfang,bende,bsw,bmonitor;   /* extern */

/* --------------------------------------------------------------- 
 Liest einen Datenblock im B2 HASYLAB Format 
 Titelzeilen sind schon weg. zuerst Kopfzeile des Blocks, am
 Ende dann 0000.0 <RET> 
--------------------------------------------------------------- */

long int blockles(FILE *in, float *px)    /* px zeigt auf x[0] */

{ int ir,schluss;
  long int n,n1,i;
  char buffer[250];
  char *form1 = "%f %f %f %f %d",
       *form2 = "%f %f %f %f";
  float temp, wert, mon, zeit, winkel;
  int blocknr;

  banfang=bsw=bende=0; 

/* erste Zeile enthaelt "b " und 5 Worte   */

  fgets(buffer,250,in);
  if ( feof(in) ) return 0 ;
  if ( strncmp(buffer,"b temp",6)!=0 ) return -1 ;

/* zweite Zeile enthaelt Blockdaten  */

  fgets(buffer,250,in);
  ir = sscanf(buffer,form1,&temp,&banfang,&bsw,
              &bende,&blocknr); 
  if ( (banfang > bende) || (banfang<=0) || (bsw==0) || (ir<5))
    {fprintf(stderr,"Blockzeile seltsam! \n >>%s<< ",buffer);
     fprintf(stderr,"Von %.3f mit %.3f bis %.3f \n",banfang,bsw,bende);
     return -1;}
  n = ((bende-banfang)/bsw +.00001) + 1;
/*  if (redeviel)
    {
    fprintf(stdout,"Von %.3f mit %.3f bis %.3f, n=%d \n",banfang,bsw,bende,n);
    } */
  if (n>MAXWERTE) {printf("Zuviele Messpunkte in Block!!!\n");
                   return -1; }

/* n Messwerte einlesen */

  n=0; schluss=0;
  while (!schluss && !feof(in))
  {
   fgets(buffer,250,in);
   if (0==strncmp(buffer,"      0.0000",12) )
      schluss=1; /* Block zuende, eventuell verfr"uht  */
   else
     {
     ir = sscanf(buffer,form2,&wert,&mon,&zeit,&winkel);
     if (  feof(in) ||  ir<4 )
       { printf("\n Fehler beim einlesen in <blockles> in Zeile ");
	 printf("%d\n",i);
	 printf("%s",buffer);
	 schluss=1; }
     else
       { *++px = wert;
	 n++; }
     } /* else */
  }
  bende=banfang+(n-1)*bsw; bmonitor=mon;
 /*       printf("-->%.3f sw=%.3f %.3f %.3f n=%d\n ",   
              banfang,bsw,bende,mon,n);  */
  return n;
}
/* --------------------------------------------------------------- */
/*  Schreiben eines Blockes in ein DUO-Datei                       */
/* --------------------------------------------------------------- */

void blockschreib(FILE *out, float y[], float start, float sw, float ziel,         		  long int n)
{ long int i;
  if (n>0)
  {
   /* Schreiben der Blockbefehlszeile fuer das DUO-File */
   fprintf(out,"B\n       0.000  %10.5f  %9.5f   %10.5f    %8.4f  %2d   %4d\n"	          ,start,sw,ziel,300.0,1,n);
   /* Schreiben der Daten */
   for (i=1; i<=n; i++)
    { fprintf(out,"%7.0f ",y[i]);
      if ((i % 8)  <= 0)
        fprintf(out,"\n");
    }
   fprintf(out,"\n");
  }
}

/* --------------------------------------------------------------- */
/*  Hauptprozedur 
    Erzeugt aus einem HASYLAB-Datei mit mehreren Bloecken 
    eine DAT-Datei, wobei die Luecken nicht aufgefuellt werden.  
    Ausserdem wird eine Protokolldatei geschrieben (.log).
    Ein vom Benutzer eingegebener Referenzmessbereich, der immer wieder
    zwischen den 'normalen' Winkelbereichen gemessen wird, wird mit seiner
    Intensitaetssumme benutzt, um im nachhinein eine interne Skalierung
    aller anderen Messbloecke vorzunehmen.
    */ 
/* --------------------------------------------------------------- */

void konvertiere(char *quelle,char *z1datei,char *z2datei)
{ FILE *in, *out1, *out2, *korr;
  fpos_t filepos;
  char buffer[250];
  float anfang, sw, ende, a, b, fdiff;
  float bwerte[MAXBLOCK][MAXWERTE], banf[MAXBLOCK], bschluss[MAXBLOCK],
        bschritt[MAXBLOCK], bmon[MAXBLOCK], x[MAXWERTE], y;
  float korranf,korrsw,korrende,korrmonitor,korrsumref,ksum;
  float korrfaktor[MAXBLOCK],k1,k2,fksum,korrsum[MAXBLOCK];
  float excvon[MAXBLOCK],excbis[MAXBLOCK];
  int nexc;
  int bmarkiert[MAXBLOCK], istkorr[MAXBLOCK];
  long int bnn[MAXBLOCK];
  long int nges, n=0, nn, i, j, k, iret, idiff;
  int ib, nblock, iberster, ibletzter, 
      ibkorr, nkorrblock, ikvor,iknach,ikorrerster;

  /*  initialisiere markierungsfeld fuer zu ignorierende Bloecke */  
  for (i=1; i<=MAXBLOCK; i++) bmarkiert[i]=0;

  printf("\nBitte gebe  ANFANG, SCHRITTWEITE, ENDE und MONITORZAHL des",
       "Referenzmess-");
  printf("\nbereiches ein, der mit seiner Intensitaetssumme zur internen ");
  printf("\nSkalierung der Messdaten verwendet werden soll:");
  printf("\n(Gebe unsinnige Werte ein, wenn keine Referenzen gemessen wurden");
  printf("\nund die Daten einfach aneinandergehaengt werden sollen)\n");
  scanf("%f %f %f %f",&korranf,&korrsw,&korrende,&korrmonitor);
/* korranf=50.25; korrsw=0.006; korrende=50.7; korrmonitor=10000; */
  printf("Referenzbereich: Von %.3f mit sw %.3f bis %.3f, Monitor=%.0f\n\n",
          korranf,korrsw,korrende,korrmonitor);
  
  if ( ((in=fopen(quelle,"rt")) == NULL) || feof(in) )
   {  fprintf(stderr, "Cannot open input file %s.\n",quelle);
      exit(1); }
  if ( (out1=fopen(z1datei,"wt")) == NULL )
   {  fprintf(stderr, "Fehler beim Oeffnen der Datei %s.\n",z1datei); 
      exit(1); }
  if ( (out2=fopen(z2datei,"wt")) == NULL )
   {  fprintf(stderr, "Fehler beim Oeffnen der Datei %s.\n",z2datei); 
      exit(1); }

/* Umkopieren der Titelzeilen */

  fprintf(out1,"Datei %s konvertiert mit zusammenskal.c\n",quelle);
  fprintf(out2,"Datei %s konvertiert mit zusammenskal.c\n",quelle);
  fprintf(out1,
     "Referenzbereich fuer interne Skalierung ist %.3f %.3f %.3f %.0f",
     korranf,korrsw,korrende,korrmonitor);
  fprintf(out2,
     "Referenzbereich: Von %.3f mit sw %.3f bis %.3f, Monitor=%.0f\n\n",
     korranf,korrsw,korrende,korrmonitor);
  i=0;
  while (!feof(in) && (i++<40))
   {  fgets(buffer,250,in);
      if (0==strncmp(buffer,"titelende",9))
         i=42;
      else 
       { fprintf(out1,"%s",buffer);fprintf(out2,"%s",buffer);} 
   }
  fprintf(out1,"\n"); 
  fprintf(out2,"\n"); 

/* ================== Hauptteil =====================================   */

/* -------- Einlesen aller Daten-Bloecke  --------------------------- */    

  /* anfang, ende,nges      gehoeren zum gesamten Datenblock,
     banfang,bsw,bende      gehoeren zum "neuen" Block      */

  nges=0; anfang=ende=0;
  ib=0; ibkorr=0; 

  while (!feof(in))
     {
     /* Einlesen des Datenblocks, Ergebnis skaliert in bwerte[]  */
     iret = blockles(in,x); 
     if (iret<0)
       { printf(
          "\n Schlimmer Fehler beim Einlesen der Daten in blockles \n"); 
         fprintf(out2,
          "\n Schlimmer Fehler beim Einlesen der Daten in blockles \n"); 
         if (feof(in))
           printf("End of File gefunden\n");
         break;
       }
     if ((iret<5)) continue; 

     /* nun hat der Block also >0 Werte */
     ib++;

     if (nges > MAXGESWERTE)
       {printf("Maximale Messpunktezahl ueberschritten -- Abbruch.");
        exit(1);}
     if ( (ib >= MAXBLOCK) ) 
       {printf("Zuviele Bloecke -- Abbruch. ");
        exit(1);}

     /* Ist der aktuelle Block ein Referenzwinkelbereich fuer die
        Intensitaetskorrektur ? */

     if (   (fabs(banfang-korranf)<0.0005) 
         && (fabs(bsw-korrsw)<0.0002)
         && (fabs(bende-korrende)<0.0005) 
         && (fabs((bmonitor-korrmonitor)/bmonitor) < 0.0001) )
          { ksum=0; for (i=1; i<=iret; i++) ksum += x[i]; 
            if (ksum>0) { istkorr[ib]=1; korrsum[ib]=ksum; ibkorr++; 
                          if (ibkorr==1) {korrsumref=ksum; ikorrerster=ib;} }
          }
     else { istkorr[ib]=0; }

     /* Eingelesenen Block merken */
     n=iret; nges+=n;
     bnn[ib]=n; bschritt[ib]=bsw; bmon[ib]=bmonitor;
     banf[ib]=banfang; bschluss[ib]=bende;
     for (i=1; i<=n; i++) bwerte[ib][i]=x[i];

     printf("Block %d: %.3f \- %.3f (sw=%.3f, Monitor=%.0f, n=%d)",
             ib,banfang,bende,bsw,bmonitor,n);
     fprintf(out2,
            "Block %d: %.3f \- %.3f (sw=%.3f, Monitor=%.0f, n=%d)",
             ib,banfang,bende,bsw,bmonitor,n);
     if ((istkorr[ib]==1)) {
       printf(      " Referenzsumme=%.0f",ksum); 
       fprintf(out2," Referenzsumme=%.0f",ksum); }
     printf("\n"); fprintf(out2,"\n");

     if ((banfang<anfang)||(ib==1)) anfang=banfang;
     if ((bende>ende)||(ib==1)) ende=bende;

    } /* Ende der While-Schleife ueber Eingabe-Datei */

     fclose(in);

/* ----------------------------------------------------- */
/*    Einlesen der Daten beendet                         */
/* ----------------------------------------------------- */

     nblock=ib; nkorrblock=ibkorr;
     printf("\n%d Referenzbloecke gefunden.",nkorrblock);
     if (nkorrblock<=1) 
        printf("\nEs wird keine interne Intensitaetseichung durchgefuehrt!\n"); 

/*------ nun enthalten 

     banf[], bschluss[]: Blockanfang 2theta, Blockende,
     bwerte[ib][ix]:     Messwerte
     bschritt[]:         Schrittweite,
     bnn[]:              Anzahl der Messwerte in diesem Block
     bmon[] :            Monitorcounts pro Schritt.
     anfang, ende:       kleinster und groesster Winkel der ganzen Datei
     nges:		 gesamte Zahl an Messpunkten in allen Bloecken.
     nblock:             gesamte Zahl von Messbloecken.
     nkorrblock:         Anzahl von gefundenen Referenzbloecken.
     istkorr[]:          0/1 je nachdem, ob dies ein zu Vergleichszwecken
                         gemessener Referenzblock ist (ja=1/nein=0).
     korrsum[]:          Intensitaetssumme im Referenzblock.
     korrsumref:         Intensitaetssumme des ersten Referenzblocks.
     ikorrerster:        Dessen Nummer.

     Man muss nun die Intensitaetssummen der Referenzbloecke in Korrektur-
     faktoren fuer die 'normalen' Messbloecke umrechnen.
     Danach koennen mehrfach gemessene Bloecke gemittelt werden, mit 
     entsprechend kleineren sigmas.
*/

/* ----------------------------------------------------------------------
        bestimme Korrekturfaktoren fuer 'normale' Messbloecke 
   --------------------------------------------------------------------*/

     for (ib=1; ib<=nblock; ib++) korrfaktor[ib]=1;  /* Inititalisierung */

     if ((nkorrblock>1) && (nkorrblock<nblock))
     { 
       printf("\nErster Referenzblock %d hat Intensitaetssumme %.0f \n",
          ikorrerster,korrsumref);
       fprintf(out2,"\nErster Referenzblock %d hat Intensitaetssumme %.0f\n",
          ikorrerster,korrsumref);

       for (ib=ikorrerster+1; ib<=nblock; ib++) 
       { /* Schleife ueber alle zu korrigierenden Messbloecke */
         if (istkorr[ib]==1) 
           {korrfaktor[ib]=korrsumref/korrsum[ib];}
         else
         {         
         /* suche vorhergehenden Referenzblock */
         ikvor=0; i=ib-1;
         while ((ikvor==0) && (i>=ikorrerster))
            { if (istkorr[i]==1) ikvor=i; i--; }
         /* suche nachfolgenden Referenzblock */
         iknach=0; i=ib+1; 
         while ((iknach==0) && (i<=nblock))
            { if (istkorr[i]==1) iknach=i; i++; }
         /* berechne Korrekturfaktor, linear interpoliert */
         if (iknach==0) korrfaktor[ib]=korrsumref/korrsum[ikvor];
         if (iknach>0)  /* interpoliere linear */
           { k1=korrsum[ikvor]; k2=korrsum[iknach];
             fksum=k1+1.0*(ib-ikvor)*(k2-k1)/(1.*(iknach-ikvor));
             korrfaktor[ib]=korrsumref/fksum;
           }
         } /* ende if ist Referenzblock */
       } /* ende for */

/* --- Korrektur der gemessenen Intensitaetswerte  ---- */

       for (ib=1; ib<=nblock; ib++)
         { printf("Korrfaktor fuer Block %3d ist %.4f", ib,korrfaktor[ib]);
           if (istkorr[ib]==1) printf(" (Referenzblock)");
           printf("\n");
           fprintf(out2,"Korrekturfaktor fuer Block %3d ist %.4f"
               ,ib,korrfaktor[ib]);
           if (istkorr[ib]==1) fprintf(out2," (Referenzblock)");
           fprintf(out2,"\n");
           if (istkorr[ib]==0)
             for (i=1; i<=bnn[ib]; i++)
               bwerte[ib][i]=bwerte[ib][i]*korrfaktor[ib];
         }
     } /* ende if Korrektur moeglich */
   
/* ---------------------------------------------------------------------
    MERGEN DER DATEN: suche zusammenhaengenden Bereich in 2 theta 
   ------------------------------------------------------------------ */

     printf("\nGesamter Datensatz: Anfang=%.3f Ende=%.3f Messpunkte=%d\n",
            anfang,ende,nges);
     fprintf(out2,
            "Gesamter Datensatz: Anfang=%.3f Ende=%.3f Messpunkte=%d\n",
            anfang,ende,nges);
     if (nges==0) {printf("Keine Daten gefunden.\n"); exit(1);}

     i=0; nn=0; ib=1;
     /* finde ersten normalen Block */
     while ((i==0) && (ib<=nblock))
     {  if (istkorr[ib]==0) i=1;
        else ib++;
     }
     iberster=ib; banfang=banf[iberster]; bende=bschluss[iberster]; 
     nn+=bnn[iberster];
     i=0; ib++; ibletzter=iberster;

     while ((i==0) && (ib<=nblock))
     {  if (istkorr[ib]==0)
          {if (banf[ib]+bschritt[ib] > bende) 
             { /* ok:  hoeherer Winkel */
             nn+=bnn[ib]; ibletzter=ib;
             bende=bschluss[ib]; 
             }
           else {i=-1;}  /* nicht ok: hier abbrechen */ 
          }
        ib++;
     }
     ib--;
     printf("Gefundener aufeinanderfolgender 2theta-Bereich:\n");
     printf("Anfang=%.3f  Ende=%.3f   Messpunkte=%6d \n",
            banfang,bende,nn);
     printf("Erster  Block: %4d \n",iberster);
     printf("Letzter Block: %4d \n",ibletzter);
     fprintf(out2,"Gefundener aufeinanderfolgender 2theta-Bereich:\n");
     fprintf(out2,"Anfang=%.3f  Ende=%.3f   Messpunkte=%6d \n",
            banfang,bende,nn);
     fprintf(out2,"Erster  Block: %4d \n",iberster);
     fprintf(out2,"Letzter Block: %4d \n",ibletzter);

/* -------- merge mehrfache Bloecke ------------------------------------ */    
     if (ibletzter < nblock ) 
     { printf(
       "Die folgenden Bloecke werden auf Doppelte im Messbereich abgesucht\n"); 
       fprintf(out2,
       "Die folgenden Bloecke werden auf Doppelte im Messbereich abgesucht\n"); 
       for (i=(ibletzter+1); i<=nblock; i++)   
       {  /* Schleife ueber restliche Bloecke */
         if (istkorr[i]==0)  
         {  /* normaler Block, keine Referenzmessung  */
         printf(" %3d ",i);
         for (j=iberster; j<=ibletzter; j++) /* Schleife ueber erste Bloecke */
         { if ( (banf[i]==banf[j])&&(bschluss[i]==bschluss[j])
             &&(bschritt[i]==bschritt[j]) && (istkorr[j]==0) )
           { /* dieser Block ist gleich! Also dazuaddieren bzw. mergen */
             printf("Addiere Block %3d zu %3d\n", i, j);
             fprintf(out2,"Addiere Block %3d zu %3d\n", i, j);
             if ( (bmon[i]<=0) || (bmon[j]<=0) )
               {printf("Monitorcounts=0!!! -- Abbruch."); exit(1);}
             for (k=1; k<=bnn[j]; k++)
                bwerte[j][k] += (bwerte[i][k]*bmon[j]/bmon[i]) ;         
             ++bmarkiert[j];  /* flag setzen, dass hier dazuaddiert wurde */ 
             break; /* aus innnerer for schleife */
            } /* end if gleich */
         } /* end for j */
         } /* end if istkorr==0 */
       }  /* end for i */
     }   /* end if */

/* ---------------------------------------------------------------------- */
/* -------- 3. Schreibe Bloecke fortlaufend in DAT Datei ---------------- */    
/* ---------------------------------------------------------------------- */

   fprintf(out1,"DAT  10 %.4f %.4f %.4f   290  1 %5d\n",
                banfang,0.002,bende,nn);
   fprintf(out2,"EXCLUDED REGIONS\n");
   printf(      "\nBloecke in Datei schreiben ... \n");
   fprintf(out2,"\nBloecke in Datei schreiben ... \n");
 
   k1=0; k2=0; nexc=0;

   for (i=iberster; i<=ibletzter; i++)
   { /* Schleife ueber alle Bloecke */
     if (istkorr[i]==0) 
     {
       printf(" %3d ",i); fprintf(out2," %3d ",i);

       if (bmarkiert[i]==0)  /* einfach gemessener Block */
         {
         for (j=1; j<=bnn[i]; j++)
           { k1=banf[i]+bschritt[i]*(j-1);
             if (!((j==1) && (k1<=k2)))
             { fprintf(out1,"%.3f %.1f %.3f %.1f \n", k1,
                    bwerte[i][j], sqrt(bwerte[i][j]), 0);
               k2=k1; /* merke letzten 2-theta Wert */
             }
           }
         }
       else  /* mehrfach gemessener Block */
         {for (j=1; j<=bnn[i]; j++)
           { k1=banf[i]+bschritt[i]*(j-1);
             if (!((j==1) && (k1==k2)))
             { n = bmarkiert[i]+1; y = bwerte[i][j]/n;
               fprintf(out1,"%.3f %.1f %.3f %.1f \n", k1, y, sqrt(y/n), 0);
               k2=k1;
             }
           }
         } /* end else */  

      /* excluded region (=fehlender Messbereich) merken, wenn schon ein 
         Block vorbei  */ 
      if (i>iberster) 
         if ( banf[i]-bende+bschritt[i] > 0) 
            nexc++; excvon[nexc]=bende+0.001;
                    excbis[nexc]=banf[i]-0.001;
      bende=bschluss[i]; /* altes Blockende merken */

      } /* end if dieser Block ist normaler Messblock */

     } /* ende schleife alle bloecke */

   /* Noch excluded regions in Logbuch schreiben */
   fprintf(out2,"\n\n");
   for (i=1; i<=nexc; i++) 
      fprintf(out2,"EXC> %8.3f %8.3f\n",excvon[i],excbis[i]);

   fclose(out1); fclose(out2);

} ; /* Ende konvertiere */


/* -------------------------------------------------------------------- */
/* >>>>>>>>>>>>>>>>>>   H A U P T P R O G R A M M   <<<<<<<<<<<<<<<<<<< */
/* -------------------------------------------------------------------- */

void main(int argc, char *argv[])
{  
   char quelle[22],z1datei[22],z2datei[22];
   int i;

   if (argc != 2)
     { printf("Bitte einen Dateinamen beim Aufruf angeben.\n"); 
       exit(0); }
   else
   {
     strcpy(quelle,argv[1]);
     printf("Quelle %s ---> ",quelle);
     /* nun Namen bilden mit neuer Endung .duo und .log */
     i=0;
     while ((quelle[i] != '.') && (quelle[i] != '\0'))
     { z1datei[i] = z2datei[i] = quelle[i];
       i++;
     }
     z1datei[i] = z2datei[i] = '.';
     z1datei[i+1] = 'z'; 
     z1datei[i+2] = 'u'; 
     z1datei[i+3] = 's'; 
     z2datei[i+1] = 'l';
     z2datei[i+2] = 'o';
     z2datei[i+3] = 'g';
     z1datei[i+4] = z2datei[i+4] = '\0';

     printf("Ziel: %s \n",z1datei);
     printf("Logbuch: %s \n",z2datei);

     konvertiere(quelle,z1datei,z2datei);  
     printf("\n FERTIG. \n");
   }
}