Thread: New to C please help

  1. #1
    Registered User
    Join Date
    Sep 2008
    Posts
    15

    New to C please help

    Hi everyone Im just starting to learn c becasue I want to edit software that I use a lot. Anyway the software is designed to analyse data for angular distributions. Its a Radware Source code called legft. It only allows for 10 angles to be fitted at any time, but I want it to fit measurements from 16 different angles. Also the out put for my measured data and the fitted data(These are .out and .prt files) is only to 2 decimal places and I need it to be to at least 4 decimal places. Anyway I am new to C and I have been trying to learn how to change the source code to enable me to do what I want.
    I have pasted the code below any help would and I wouldbe extremly gratefull. Cheers
    Code:
    #include <stdio.h>
    #include <stdlib.h>
    #include <unistd.h>
    #include <math.h>
    #include <string.h>
    
    #include "util.h"
    
    int main(void)
    {
      double d1, d2;
      double enc, enc2, enc4, enc6, enc8;
      float  r1, r2, a[3][7], b[10][17];
      float  e[17], g[7][7], p[10][17], r, angle[17], w[17], y[3][17];
      float  thick[17], trans, wnorm[17], x2[3], fk[7];
      float  egamma, pq[7], ewnorm[17], deg, err[3][7], att;
      int    iout1, i, j, jj, nt, nr, ir, irr, nnt, nrs, nrt;
      char   nf[80], of2[80] = "legft.prt", of3[80] = "legft.out";
      char   line[120], *nd;
      FILE   *file1, *file2, *file3;
      /* ************************************
         program legft
         program to fit angular distributions with legendre polynomials
         by the least squares method
         author:           D.C. Radford
         first written:    Feb 1989
         latest revision:  Jan 8 1992
         description:
         This program fits angular distributions using legendre polynomials
         It does two fits to each data set; the first fit is with two
         (P0+P2), the second is with three terms (P0+P2+P4). The
         input file is a modified .eng file, created by running the program
         ENERGY and editing its output file.
         The 1st line must contain the number of angles
         The 2nd line must contain the angles in degrees
         The 3rd line must contain normalisation values
         The 4th line must contain the errors on these values
         The fifth line is used for an optional correction for gamma-ray
         absorbtion in the target backing. If you do not want this correction,
         leave a blank line. Otherwise put the target backing thickness seen
         at each angle (in arbitrary units). At the end of
         each data set, at the end of the line beginning EGAMMA =, you
         must then add the absorbtion per unit thickness for that gamma ray
         (Format F10.5).
         Output is to 2 files; legft.prt has a brief summary of the results
         and is printed automatically, and legft.out has more detail and can
         be printed by the user.
         **************************************** */
    
      printf(
        "\n\n"
        "Welcome to legft.\n\n"
        "This program fits angular distributions using legendre polynomials.\n"
        "It does two fits to each data set;\n"
        " the first fit is with two terms (P0+P2),\n"
        " the second is with three terms (P0+P2+P4).\n"
        "The expected input file is a modified .eng file, created by running\n"
        " the program ENERGY and editing its output file.\n"
        "The first line must contain the number of angles        (INTEGER)\n"
        "The second line must contain the angles in degrees      (up to 10 REALs)\n"
        "The third line must contain normalisation values        (up to 10 REALs)\n"
        "The fourth line must contain the errors on these values (up to 10 REALs)\n"
        "The fifth line is used for an optional correction for\n"
        " gamma-ray absorbtion in the target backing.  If you do not want\n"
        " this correction, leave a blank line.  Otherwise put the target\n"
        " backing thickness seen at each angle (in arbitrary units,\n"
        " up to 10 REALs). At the end of each data set, at the end of the line\n"
        " beginning EGAMMA =...., you must then add the absorbtion per unit\n"
        " thickness for that gamma ray      (REAL)\n\n"
        "        D. C. Radford    Feb. 1989.\n");
      file2 = open_pr_file(of2);
      file3 = open_pr_file(of3);
    
      /* ask for data file name */
     START:
      while (1) {
        cask("Input file name = ? (default .ext = .eng)", nf, 80);
        setext(nf, ".eng", 80);
        if ((file1 = open_readonly(nf))) break;
      }
      fgets(line, 120, file1);
      nt = strtol(line, &nd, 0);
      printf("%d angles...\n", nt);
      if (nt < 1 || nt > 17) {
        printf("No. of angles must be > 0 and < 18\n");
        goto NEXTFIL;
      }
      for (i = 0; i < nt; i++) {
        angle[i] = ewnorm[i] = thick[i] = 0.0;
        wnorm[i] = 1.0;
      }
    
    #define R(data) fgets(line, 120, file1); nd = line; if (strlen(line)>=2*nt) { \
                    for (i = 0; i < nt; i++) {data[i] = strtod(nd, &nd); nd++;}}
      R(angle);
      R(wnorm);   /* wnorm = normalisation intensity */
      R(ewnorm);  /* ewnorm = normalisation error */
      R(thick);   /* thick = thickness of target backing */
    #undef R
    
      fprintf(file2,
    	  "\f\n File = %s\n\n"
    	  "   Egamma +- error          a0 +- error    a2/a0 +- error"
    	  "    a4/a0 +- error     Chisq\n", nf);
      iout1 = 0;
     NEXTDAT:
      for (nnt = 0; nnt < nt; ++nnt) {
        fgets(line, 120, file1);
        if (sscanf(line, "%*d%*f%*f%f%f", &w[nnt], &e[nnt]) != 2) goto NEXTFIL;
        /* w,e = area,error */
        if (w[nnt] == 0.0f && e[nnt] == 0.0f) nnt--;
      }
      fgets(line, 120, file1);
      att = 0.0f;
      if (sscanf(line, "%*14c%f%*3c%f%*21c%f", &egamma, &deg, &att) < 2)
        goto NEXTFIL;
      fgets(line, 120, file1);
      
      /* att = attenuation of gamma per unit thickness (c.f. thick) */
      trans = 1.0f - att;
      for (nnt = 0; nnt < nt; ++nnt) {
        r1 = e[nnt] / w[nnt];
        r2 = ewnorm[nnt] / wnorm[nnt];
        e[nnt] = sqrt(r1 * r1 + r2 * r2);
        d1 = (double) trans;
        d2 = (double) thick[nnt];
        w[nnt] /= wnorm[nnt] * pow(d1, d2);
        e[nnt] *= w[nnt];
      }
      /* calculate legendre polynomials */
      for (nnt = 0; nnt < nt; ++nnt) {
        enc  = cos(angle[nnt] * 0.0174533f);
        enc2 = enc*enc;
        enc4 = enc2*enc2;
        enc6 = enc4*enc2;
        enc8 = enc4*enc4;
        p[nnt][0] = 1.0f;
        p[nnt][2] = (enc2 * 3.0 - 1.0) / 2.0;
        p[nnt][4] = (enc4 * 35.0 - enc2 * 30.0 + 3.0) / 8.0;
        p[nnt][6] = (enc6 * 231.0 - enc4 * 315.0f + enc2 * 105.0 - 5.0) / 16.0;
        p[nnt][8] = (enc8 * 6435.0 - enc6 * 12012.0 + enc4 * 6930.0 -
    		 enc2 * 1260.0 + 35.0) / 128.0;
      }
    
      for (nr = 3; nr <= 5; nr += 2) {
        jj = nr/2 - 1;
        for (irr = 0; irr < nr; irr += 2) {
          fk[irr] = 0.f;
          for (ir = 0; ir < nr; ir += 2) {
    	b[irr][ir] = 0.f;
          }
        }
        for (nnt = 0; nnt < nt; ++nnt) {
          for (irr = 0; irr < nr; irr += 2) {
    	for (ir = 0; ir < nr; ir += 2) {
    	  b[irr][ir] += p[nnt][ir]*p[nnt][irr] / (e[nnt]*e[nnt]);
    	}
    	fk[irr] += w[nnt] * p[nnt][irr] / (e[nnt]*e[nnt]);
          }
        }
        g[0][0] = 1.f / b[0][0];
        for (ir = 2; ir < nr; ir += 2) {
          r = b[ir][ir];
          for (i = 0; i < ir; i += 2) {
    	pq[i] = 0.f;
    	for (j = 0; j < ir; j += 2) {
    	  pq[i] += g[j][i] * b[ir][j];
    	}
    	r -= pq[i] * b[ir][i];
          }
          g[ir][ir] = 1.f / r;
          for (i = 0; i < ir; i += 2) {
    	g[ir][i] = g[i][ir] = -pq[i] / r;
    	for (j = 0; j < ir; j += 2) {
    	  g[j][i] += pq[i] * pq[j] / r;
    	}
          }
        }
        for (nrt = 0; nrt < nr; nrt += 2) {
          err[jj][nrt] = sqrt(g[nrt][nrt]);
          a[jj][nrt] = 0.f;
          for (nrs = 0; nrs < nr; nrs += 2) {
    	a[jj][nrt] += g[nrs][nrt] * fk[nrs];
          }
        }
        for (nnt = 0; nnt < nt; ++nnt) {
          y[jj][nnt] = 0.f;
        }
        for (nnt = 0; nnt < nt; ++nnt) {
          for (nrt = 0; nrt < nr; nrt += 2) {
    	y[jj][nnt] += a[jj][nrt] * p[nnt][nrt];
          }
        }
        x2[jj] = 0.f;
        for (nnt = 0; nnt < nt; ++nnt) {
          r1 = (y[jj][nnt] - w[nnt]) / e[nnt];
          x2[jj] += r1 * r1;
        }
        x2[jj] /= (float) (nt - (nr + 1) / 2);
        for (nrt = 2; nrt < nr; nrt += 2) {
          a[jj][nrt] /= a[jj][0];
          err[jj][nrt] /= a[jj][0];
          r1 = err[jj][1] * a[jj][nrt] / a[jj][0];
          err[jj][nrt] = sqrt(err[jj][nrt]*err[jj][nrt] + r1*r1);
        }
      }
      fprintf(file3, "    Angle  cos**2         Data +- error"
    	  "      Fit(p0+p2)   Fit(p0+p2+p4)\n");
      for (nnt = 0; nnt < nt; ++nnt) {
        enc = cos(angle[nnt] * .0174533f);
        enc2 = enc*enc;
        fprintf(file3, "%8.2f %7.4f %12.2f %8.2f %15.2f %15.2f\n",
        angle[nnt], enc2, w[nnt], e[nnt], y[0][nnt], y[1][nnt]);
      }
      fprintf(file3, "   Chisq = %8.3f %8.3f\n", x2[0], x2[1]);
      fprintf(file3,
    	  "    File = %s  a0 = %.2f +- %.2f, [%.2f +- %.2f]\n"
    	  "  Egamma = %.4f +- %.4f   a2/a0 = %.4f +- %.4f, [%.4f +- %.4f]\n"
    	  "  [a4/a0 = %.4f +- %.4f]\n\n",
    	  nf, a[0][0], err[0][0], a[1][0], err[1][0], egamma, deg,
    	      a[0][2], err[0][2], a[1][2], err[1][2], a[1][4], err[1][4]);
      fprintf(file2,
    	  "%10.3f %7.3f %12.2f %7.2f %9.3f %7.3f %27.3f\n"
    	  "%31.2f %7.2f %9.3f %7.3f %9.3f %7.3f %9.3f\n",
    	  egamma, deg,
    	  a[0][0], err[0][0], a[0][2], err[0][2], x2[0],
    	  a[1][0], err[1][0], a[1][2], err[1][2], a[1][4], err[1][4], x2[1]);
      ++iout1;
      if (iout1 >= 28) {
        printf("\f\n File = %s\n\n"
    	   "   Egamma +- error          a0 +- error    a2/a0 +- error"
    	   "    a4/a0 +- error     Chisq\n", nf);
        iout1 = 0;
      }
      goto NEXTDAT;
    
     NEXTFIL:
      fclose(file1);
      if (caskyn("Would you like to process more files? (Y/N)")) goto START;
      pr_and_del_file(of2);
      fclose(file3);
      printf("Condensed output sent to printer.\n"
    	 "For more detail, print file legft.out.\n");
      return 0;
    } /* legft */

  2. #2
    Kernel hacker
    Join Date
    Jul 2007
    Location
    Farncombe, Surrey, England
    Posts
    15,677
    This is what I would say to the original author (I'm about 16 years too late tho'):
    Let me guess: You used to program in Basic? That's why you so happily use goto, when do-while or while would be perfectly fine choices.

    Code:
    #define R(data) fgets(line, 120, file1); nd = line; if (strlen(line)>=2*nt) { \
                    for (i = 0; i < nt; i++) {data[i] = strtod(nd, &nd); nd++;}}
      R(angle);
      R(wnorm);   /* wnorm = normalisation intensity */
      R(ewnorm);  /* ewnorm = normalisation error */
      R(thick);   /* thick = thickness of target backing */
    #undef R
    You should really do this as a function. There is ABSOLUTELY no need to use a multi-line macro here.

    Code:
      float  e[17], g[7][7], p[10][17], r, angle[17], w[17], y[3][17];
    Somewhat longer names would help. (angle is fine).
    .

    You may want to read up on printf format strings. Also, if you want more precise results, you probably need to use double instead of float. This will affect both the variable declarations and the scanf() that you read the data with.

    --
    Mats
    Compilers can produce warnings - make the compiler programmers happy: Use them!
    Please don't PM me for help - and no, I don't do help over instant messengers.

Popular pages Recent additions subscribe to a feed