home *** CD-ROM | disk | FTP | other *** search
/ Meeting Pearls 1 / Meeting Pearls Vol 1 (1994).iso / installed_progs / gfx / lise2.1 / lise / src / rvt2.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-03-31  |  11.6 KB  |  409 lines
  1. /*
  2.   split one spectra file, containing 4 or 8 spectra, into
  3.   single spectra, possibly reverse spectra, subtract Background
  4.   Calculate PAC ratio function:
  5.  
  6.                                         1/4
  7.           2   (product(Si(180) - BACKi))
  8.    R(t) = - * --------------------------
  9.           3   (product(Sj(90) - BACKj))^(1/4)
  10.  
  11.   Input is read from the file specified.
  12.   The peakpositions, phase and parity are read from *.def or global.def
  13. */
  14.  
  15. #include <stdio.h>
  16. #include <spec.h>
  17. #define TMPFILE "global.def"
  18.  
  19. char deffile[80];
  20. float *spc, *err, *tim, *usp1, *usp2, *uer1, *uer2;
  21. int range=1024,
  22.     tri=0,
  23.     flg_v=FALSE,
  24.     flg_ba=0,
  25.     lastn=10,
  26.     dev=50;
  27.  
  28. int *back;
  29.  
  30. help()
  31. {
  32.   printf("Splitting spectrum file into single spectra, possibly reversing\n");
  33.   printf("them, subtracting Background\n");
  34.   printf("Calculating PAC Ratio function:\n");
  35.   printf("                                        1/4       \n");
  36.   printf("        2   (product(Si(180,t) - BACKi))          \n");
  37.   printf(" R(t) = - * --------------------------------- - 1 \n");
  38.   printf("        3   (product(Sj(90,t) - BACKj))^(1/4)     \n");
  39.   printf("\n");
  40.   printf("               1     2   [ 1                               \n"); 
  41.   printf(" Error(t) = ------ * - * [-- * sum(1/(Si(180,t)-BACKi)) +  \n");
  42.   printf("            R(t)+1   3   [16                               \n");
  43.   printf("                                                          1/2 \n");
  44.   printf("                            1                            ]    \n"); 
  45.   printf("                         + -- * sum(1/(Sj(180,t)-BACKj)) ]    \n");
  46.   printf("                           16                            ]    \n");
  47.   printf("\n");
  48.   printf("RVT2 file [-o file] [-back n] [-last n]\n");
  49.   printf("options:\n");
  50.   printf("   -o  file     output spectrum instead of stdout\n");
  51.   printf("   -back n1,n2,n3,..  background for all spectra\n");
  52.   printf("   -last n      takes arithmetic mean of last n channels as\n");
  53.   printf("                background. Default is, to find the background\n");
  54.   printf("                by fitting an exponential curve to the spectrum\n");
  55.   printf("   -v           verbose mode. Prints fitted background\n");
  56.   printf("   -def file    use spectra definition file different from %s\n",TMPFILE);
  57.   printf("The definition of the spectrum is read from %s\n",deffile);
  58.   printf("This file should be in a format as generated by the FPROMPT \n");
  59.   printf("routine and is organized as follows:\n");
  60.   printf("Each line consists of 4 integers.\n");
  61.   printf("The first specifies the Time 0 position.\n");
  62.   printf("The second can be +1 -1 +2 or -2 and specifies, if it is\n");
  63.   printf("     to be mirrored (-) and if it is a 0 degree (1)\n");
  64.   printf("     or 90 degree (2) spectrum\n");
  65.   printf("     If any other number is given, the spectrum is ignored\n");
  66.   printf("The last two numbers are used to specify the range\n");
  67.   printf("of the usable data area. Up to now, these parameters\n");
  68.   printf("are only significant in the last line !\n");
  69.   printf("\n(C) Rainer Kowallik\n");
  70.   exit(0);
  71. }
  72.  
  73. float sd(x,y)
  74. float x,y;
  75. {
  76. if(y!=0.0) return((float)(x / y));
  77. return(0.0);
  78. }
  79.  
  80. float root(n,x)
  81. float x;
  82. int n;
  83. {
  84. float y;
  85.  
  86. y = fabs(exp((double)(log(fabs((double)x))/((double)n))));
  87. return(y);
  88. }
  89.  
  90. main(argc,argv)
  91. int argc;
  92. char *argv[];
  93. {
  94. int  xbeg,n,m,i,max,peakptr,pacmax,peaks[30],phase[30];
  95. int  wurzel1,wurzel2,merk;
  96. char *z,*s,*comment,*defnam;
  97. FILE *fp;
  98.  
  99.    if(argc == 1) help();
  100.  
  101.    z = (char *) malloc(256);
  102.    s = (char *) malloc(256);
  103.    comment = (char *) malloc(256);
  104.    defnam = (char *) malloc(256);
  105.    back = (int *) calloc(32,sizeof(int));
  106.  
  107.    strcpy(deffile,TMPFILE);
  108.    flg_ba=0;            /* defaults to exponential fit */
  109.    for(n = 0; n < 32; n++) back[n] = -1;
  110.    flg_v=FALSE;
  111.  
  112.    if(checkopt(argc,argv,"-def",z)) strcpy(deffile,z);
  113.    strcpy(defnam,deffile);
  114.    if(checkopt(argc,argv,"-back",z)) {
  115.       i = 0; n = 0; m = 0; flg_ba=1;
  116.       while(z[n] != 0) {
  117.          s[m++] = z[n++];
  118.          if(z[n] == ',') {
  119.             s[m] = 0; m = 0; n = n + 1;
  120.             back[i++] = atoi(s);
  121.          }
  122.       }
  123.       while(i < 32) back[i++] = atoi(s);
  124.    }
  125.    if(checkopt(argc,argv,"-last",z)) {
  126.       flg_ba=2; lastn = atoi(z);
  127.    }
  128.    if(checkopt(argc,argv,"-v",z)) flg_v=TRUE;
  129.  
  130.  
  131.    spc=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  132.    err=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  133.    tim=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  134.    usp1=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  135.    usp2=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  136.    uer1=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  137.    uer2=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  138.  
  139. /* -----------------------------------------------
  140.    reading spectrum 
  141.    ----------------------------------------------- */
  142.  
  143.    max=readspec(argv[1],spc,err,tim,comment);
  144.  
  145. /* -----------------------------------------------
  146.    generating error spectrum 
  147.    ----------------------------------------------- */
  148.  
  149.    for(n=0;n<max;n++) err[n] = sqrt(spc[n]);
  150.  
  151. /* -----------------------------------------------
  152.    reading spectra definition file
  153.    ----------------------------------------------- */
  154.  
  155.    fp=fopen(defnam,"r");
  156. #ifdef UNIX
  157.    if(fp==NULL) {
  158.       home = (char *) getenv("HOME");
  159.       if(home != NULL) {strcpy(defnam,home); strcat(defnam,"/"); }
  160.       strcat(defnam,deffile);
  161.       fp=fopen(defnam,"r");
  162.    }
  163.    if(fp==NULL) {
  164.       home = (char *) getenv("LISEPRG");
  165.       if(home != NULL) {strcpy(defnam,home); strcat(defnam,"/"); }
  166.       strcat(defnam,deffile);
  167.       fp=fopen(defnam,"r");
  168.    }
  169. #endif
  170. #ifdef AMIGA
  171.    if(fp==NULL) {
  172.       strcpy(defnam,"LISE:"); strcat(defnam,deffile);
  173.       fp=fopen(defnam,"r");
  174.    }
  175. #endif
  176.    if(fp == NULL) {
  177.       printf("rvt2: could not open >%s<\n",deffile);
  178.       exit(0);
  179.    }
  180.    peakptr=0;
  181.    while(!feof(fp)) {
  182.       fscanf(fp,"%d %d %d %d\n",&peaks[peakptr],&phase[peakptr],&i,&pacmax);
  183.       peakptr=peakptr+1;
  184.    }
  185.    fclose(fp);
  186.  
  187. /* -------------------------------------------------
  188.    initializing usp1, usp2 ,uer1 and uer2
  189.    ------------------------------------------------- */
  190.  
  191.    for(n=0;n<pacmax;n++) {
  192.       usp1[n]=1.0;
  193.       usp2[n]=1.0;
  194.       uer1[n]=0.0;
  195.       uer2[n]=0.0;
  196.    }
  197.  
  198. /* ----------------------------------------------
  199.    here we start calculating the new sum spectra 
  200.   -----------------------------------------------*/
  201.  
  202.    wurzel1 = 0; wurzel2 = 0;
  203.  
  204.    for(n=0;n<peakptr;n++) {
  205.       xbeg = peaks[n];
  206.       m = phase[n];
  207. /* if negative, go and mirror spectrum */
  208.       if(m < 0) {
  209.          for(i=pacmax/2 ; i<pacmax ; i++) {
  210.             merk = spc[xbeg - i];
  211.             spc[xbeg - i] = spc[xbeg + i - pacmax];
  212.             spc[xbeg + i - pacmax] = merk;
  213.             err[xbeg - i] = sqrt(spc[xbeg - i]);
  214.             err[xbeg + i - pacmax] = sqrt(spc[xbeg + i - pacmax]);
  215.          }
  216.          xbeg = xbeg - pacmax; m = 0 - m;
  217.       }
  218.  
  219.       SubBackground(spc,err,xbeg,xbeg+pacmax,n);
  220.  
  221.       switch(m) {   /* which spectrum to add ? */
  222.       case 1:       /* 180 degree spectrum */
  223.          for(i=0;i<pacmax;i++) {
  224.          usp1[i] = usp1[i] * spc[xbeg+i];
  225.          uer1[i] = uer1[i] + sd(1.0,spc[xbeg+i]);
  226.          }
  227.          wurzel1 = wurzel1 + 1;
  228.          break;
  229.       case 2:       /* 90 degree spectrum */
  230.          for(i=0;i<pacmax;i++) {
  231.          usp2[i] = usp2[i] * spc[xbeg+i];
  232.          uer2[i] = uer2[i] + sd(1.0,spc[xbeg+i]);
  233.          }
  234.          wurzel2 = wurzel2 + 1;
  235.          break;
  236.       }
  237.    }
  238.  
  239. /* --------------------------------------------------------------
  240.          calculate ratio function and error spectrum 
  241.    -------------------------------------------------------------- */
  242.  
  243.    for(i=0;i<pacmax;i++) {
  244.       spc[i] = sd(2.0 * root(wurzel1,usp1[i]) , 3.0 * root(wurzel2,usp2[i]));
  245.       spc[i] = spc[i] - 0.666667;
  246.       err[i] = sd(1.0,1.0 + spc[i]) * (2.0 / 3.0);
  247.       err[i] = err[i] *
  248.          sqrt(sd(1.0,(float)(wurzel1 * wurzel1)) * uer1[i] + 
  249.               sd(1.0,(float)(wurzel2 * wurzel2)) * uer2[i]);
  250.    }
  251.  
  252.  
  253. /* -------------------------------------------
  254.    now we can go and save the new spectrum 
  255.    ------------------------------------------- */
  256.  
  257.    strcpy(z,comment);
  258.    n=instr("|",comment);
  259.    if(n>0) midstr(z,comment,0,n-1);
  260.    strcat(z," | RVT(PAC)");
  261.    writespec("",spc,err,pacmax,2,z);
  262.    free(uer1); free(uer2); free(usp1); free(usp2);
  263.    free(spc); free(err); free(tim);
  264.    free(back); free(defnam); free(comment);
  265.    free(s); free(z);
  266.    exit(0);
  267. }
  268. /* --------------------------------------------
  269.    finding and subtracting correct background
  270.    -------------------------------------------- */
  271. SubBackground(spc,err,xbeg,xend,j)
  272. float *spc, *err;
  273. int xbeg, xend, j;
  274. {
  275. int n,ba;
  276.  
  277.    switch(flg_ba) {
  278.    case 0:
  279.       ba = expbafit(spc,err,xbeg,xend);
  280.       break;
  281.    case 1: /* allready done when parsing arguments */
  282.       ba = back[j];
  283.       break;
  284.    case 2:
  285.       ba = sumspc(spc,xend-lastn,xend) / lastn;
  286.       break;
  287.    }
  288.    if(flg_v) printf("background=%d\n",ba);
  289.  
  290. /*    Subtracting background  */
  291.    for(n=xbeg ; n<xend ; n++) spc[n] = spc[n] - ba;
  292.    return(0);
  293. }
  294.  
  295. sumspc(spc,a,z)
  296. float spc[];
  297. int a,z;
  298. {
  299. int i,erg;
  300.  
  301.    erg=0;
  302.    for(i=a;i<z;i++) erg=erg+spc[i];
  303.    return(erg);
  304. }
  305.  
  306. /* ---------------------------------------------------------
  307.    Linear regression after logarithm of spectrum
  308.    This is only to calculate the chi**2, all
  309.    other values are lost.
  310.    UP TO NOW, THE ERRORBARS ARE NOT INCLUDED !
  311.    I HAVE FORGOTTEN MY OLD SUPERBASIC SOURCE AT HOME, SORRY.
  312.    --------------------------------------------------------- */
  313.  
  314. float linreg(spc,err,n,a,z)
  315. float spc[],err[];
  316. int n,a,z;
  317. {
  318. int i;
  319. float ysoll,x,y,sumx,sumxx,sumy,sumyy,sumxy,nges,
  320.       steigung,offset,erg,
  321.       *ylog;
  322.  
  323. ylog=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  324.  
  325. nges = (float) (z-a);    /* this we know for shure */
  326. sumx = 0.0;
  327. sumxx = 0.0;
  328. sumy = 0.0;
  329. sumyy = 0.0;
  330. sumxy = 0.0;
  331.  
  332. for(i=a;i<z;i++) {
  333.    x = (float) i;
  334.    y = spc[i];
  335.    y = y - ((float) n);
  336.    if(y<0.1) y=1.0;                 /* trap overflow on log */
  337.    y = (float) log((double) y);     /* linearization of exp function */
  338.    ylog[i] = y;                     /* store for later use */
  339.    sumx = sumx + x;
  340.    sumxx = sumxx + (x*x);
  341.    sumy = sumy + y;
  342.    sumyy = sumyy + (y*y);
  343.    sumxy = sumxy + (x*y);
  344. }
  345. sumx = sumx / nges;
  346. sumy = sumy / nges;
  347. sumxx = sumxx - (nges * sumx * sumx);
  348. sumxy = sumxy - (nges * sumx * sumy);
  349. sumyy = sumyy - (nges * sumy * sumy);
  350. steigung = sumxy / sumxx;
  351. offset = sumy - (steigung * sumx);
  352.  
  353. /* calculating chi**2 */
  354.  
  355. erg=0.0;
  356. for(i=a;i<z;i++) {
  357.    x = (float) i;
  358.    y = spc[i];
  359.    y = ylog[i];
  360.    ysoll = (steigung * x) + offset;
  361.    y = y - ysoll;
  362.    erg = erg + (y * y); 
  363. }
  364. free(ylog);
  365. return(erg);
  366. }
  367.  
  368.    
  369. /* ---------------------------------------------------------
  370.    Fitting background to exponential function.
  371.    This is done, by calculating the chi**2 with
  372.    linear regression for the logarithm of the spectrum.
  373.    We than look for a minimum of chi**2 for different
  374.    backgrounds.
  375.    UP TO NOW, THE ERRORBARS ARE NOT INCLUDED !
  376.    I HAVE FORGOTTEN MY OLD SUPERBASIC SOURCE AT HOME, SORRY.
  377.    --------------------------------------------------------- */
  378. expbafit(spc,err,a,z)
  379. float spc[],err[];
  380. int a,z;
  381. {
  382. int i,n,bamin,bamax,babest,inc,xraster,xdepth;
  383. float chisq,chibest;
  384.  
  385. xraster = 10;
  386. xdepth = 5;
  387. chibest = 1E10;
  388. bamax = sumspc(spc,z-10,z) / 10 ; /* get maximum background */
  389. bamin = bamax / 2;                /* get minimum background */
  390. for(i=1;i<xdepth;i++) {
  391.    inc = (bamax - bamin) / xraster;
  392.    if(inc<1) break;
  393.    n = bamin; while(n<=bamax) {
  394.       chisq = linreg(spc,err,n,a,z); /* just calculate chi**2 */
  395.       if(chisq < chibest) {
  396.      chibest = chisq;
  397.      babest = n;
  398.      if(flg_v) {
  399.      printf("new best: background = %d   chi**2 = %f\n",babest,chibest);
  400.      }
  401.       }
  402.       n = n + inc;
  403.    }
  404.    bamin = babest - inc;
  405.    bamax = babest + inc;
  406. }
  407. return(babest);
  408. }
  409.