• 09-11-2008
zidangus
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
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);
}
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 */

• 09-11-2008
matsp
This is what I would say to the original author (I'm about 16 years too late tho'):
Quote:

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