|
|
|
|
Programs for Theoretical Calculations
The following program is called linearize.c and it takes in a solid black image of the TMS coil under investigation and converts it into an outline. That is, it produces a single loop coil of the TMS coil. THe outline image is in units of pixels. In order to convert the units into useful units a conversion factor must be used. The factor depends on the number of dpi the coil was scanned with. In our case, the conversion factor used was 280/3.4 since the image had a diameter of 280 pixels for the diameter of 3.4cm. The call to this program is
./linearize <pgm image filename> output filename
The following is the actual code for the program.
#include <stdio.h>
#include <math.h>
#include <stdarg.h>
#include <string.h>
#include <stdlib.h>
const int STEPS = 256;
void die (char *fmt,
...)
{
va_list ap;
va_start(ap, fmt);
vfprintf (stderr,
fmt, ap);
va_end(ap);
exit (1);
}
void warn (char *fmt,
...)
{
va_list ap;
va_start(ap, fmt);
vfprintf (stderr,
fmt, ap);
va_end(ap);
}
int main (int argc, char
**argv)
{
char line [256];
int x, y;
char **image;
int i, j, k,
r;
double sum_x,
sum_y, sum;
double center_x,
center_y;
double factor_x,
factor_y;
double angle;
double radius
[STEPS];
FILE *fd;
fgets (line, 256,
stdin);
if (strncmp (line,
"P5", 2))
die
("This is not a binary PGM file ... are you trying to poison me?\n"
"read \"%s\"!\n", line);
fgets (line,
256, stdin);
if (sscanf (line,
"%d %d", &x, &y) != 2)
die
("Couldn't read the dimensions\n");
fprintf (stderr,
"The image is %d x %d pixels\n", x, y);
fgets (line,
256, stdin);
image = (char
**) malloc (y * sizeof (char *));
for (i = 0; i
< y; i++) {
image
[i] = (char *) malloc (x * sizeof (char));
}
for (i = 0; i
< y; i++)
if
(fread (image [i], 1, x, stdin) != x)
die ("Error reading the image at row %d\n", i);
if ((fd = fopen
("solid.dat", "w")) == NULL) die ("Couldn't open the solid file\n");
for (i = 0; i
< y; i++)
for
(j = 0; j < x; j++)
if (image[i][j] == 0) {
fprintf (fd, "%d %d\n", j, i);
sum_x += j;
sum_y += i;
sum++;
}
fclose (fd);
center_x = sum_x
/ sum;
center_y = sum_y
/ sum;
fprintf (stderr,
"center of mass is at (%g, %g)\n", center_x, center_y);
if ((fd = fopen
("outline.dat", "w")) == NULL) die ("Couldn't open the outline file\n");
for (k = 0; k
< STEPS; k++) {
angle
= k * 2 * M_PI / STEPS;
factor_x
= cos (angle);
factor_y
= sin (angle);
for
(r = 1; r < x; r++) {
i = center_y + r * factor_y;
j = center_x + r * factor_x;
if (image[i][j] != 0)
break;
}
printf
("%g\t%d\n", angle, r);
radius
[k] = r;
fprintf
(fd, "%d %d\n", j, i);
}
fclose (fd);
fd = popen ("gnuplot",
"w");
fprintf (fd,
"plot 'solid.dat' with dots, 'outline.dat' with line\n");
fprintf (fd,
"!xmessage Hello World!\n");
fclose (fd);
return 0;
}
The following program
is called render.c and it takes the coil image outputted by linearize.c
and puts it into the proper orientation required for the experimenter.
That is, it is capable of producing a 2 lobe system where one lobe is simply
the mirror image of the other. The following is the call to the program:
./render x1 y1 rotation1 radius1 0 x2 y2 rotation2 radius2 0 <filename of output from linearize> output name
The following is the actual code for the program.
#include <stdio.h>
#include <math.h>
#include <stdarg.h>
#include <stdlib.h>
void die (char *fmt,
...)
{
va_list ap;
va_start(ap, fmt);
vfprintf (stderr,
fmt, ap);
va_end(ap);
exit (1);
}
void warn (char *fmt,
...)
{
va_list ap;
va_start(ap, fmt);
vfprintf (stderr,
fmt, ap);
va_end(ap);
}
struct point {
double angle;
double radius;
struct point
*next;
} *list = NULL;
void add_point (double
angle, double radius)
{
struct point
*ptr;
ptr = (struct
point *) malloc (sizeof (struct point));
ptr->next = list;
ptr->angle =
angle;
ptr->radius =
radius;
list = ptr;
}
int main (int argc, char
**argv)
{
char buf[256];
double angle,
radius;
int i, invert;
double x, y,
rotation;
struct point
*ptr;
if ((argc % 5)
!= 1)
die
("Usage: %s [ x1 y1 rotation radius ] ...\n", argv[0]);
while (fgets
(buf, 256, stdin) != NULL) {
sscanf
(buf, "%lf %lf", &angle, &radius);
add_point
(angle, radius);
}
for (i = 1; i
< argc; i += 5) {
if
(sscanf (argv[i], "%lf", &x) != 1) die ("problem parsing from %d ...
x\n", i);
if
(sscanf (argv[i+1], "%lf", &y) != 1) die ("problem parsing from %d
... y\n", i);
if
(sscanf (argv[i+2], "%lf", &rotation) != 1) die ("problem parsing from
%d ...
rotation\n", i);
if
(sscanf (argv[i+3], "%lf", &radius) != 1) die ("problem parsing from
%d ...
radius\n", i);
if
(sscanf (argv[i+4], "%d", &invert) != 1) die ("problem parsing from
%d ...
radius\n", i);
for
(ptr = list; ptr; ptr = ptr->next) {
angle = (invert) ? ptr->angle : -ptr->angle;
angle += rotation;
printf ("%g\t%g\n",
x + radius * ptr->radius * cos (angle),
y + radius * ptr->radius * sin (angle));
}
printf
("\n");
}
return 0;
}
The following is the computational engine of the
program. It takes the final coil configuration output from render.c and
calculates the magnitude of the magnetic field for each point in space
above the coil at a specified distance. The current needs to be specified
within the program. The call to this program is made with the following
command:
./magnetic-field distance <output from render> datafile name
where the distance is given in pixel units determined from the conversion factor as described previously. The following is the actual code for this program.
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <math.h>
const double mu_naught
= 7.39e-3; /* T-pixels/A Dielectric */
const double current
= 50400; /* Amperes */
struct point {
double x;
double y;
};
struct segment {
struct point
start, end;
double x, y;
struct segment
*next;
} *segments = NULL;
void die (char *fmt,
...)
{
va_list ap;
va_start(ap, fmt);
vfprintf (stderr,
fmt, ap);
va_end(ap);
exit (1);
}
void warn (char *fmt,
...)
{
va_list ap;
va_start(ap, fmt);
vfprintf (stderr,
fmt, ap);
va_end(ap);
}
void add_segment (struct
point start, struct point end)
{
struct segment
*ptr = (struct segment *) malloc (sizeof (struct segment));
if (ptr == NULL)
die ("Couldn't malloc memory for segment\n");
ptr->start =
start;
ptr->end = end;
ptr->next = segments;
segments = ptr;
}
void field (double z_measure)
{
#define SQR(x) ((x)*(x))
struct segment
*ptr;
double shortest_distance;
double a_1, a_2,
a_3;
double b_1, b_2,
b_3;
double x, y,
z = z_measure;
double R2_A,
R2_B;
double L2;
double alpha1,
alpha2;
double term1,
term2, term3;
double B_x, B_y,
B_z;
double x_hat,
y_hat, z_hat, cat_in_the_hat, the_gnu_who_knew;
for (y = -200;
y < 300; y += 10) {
for
(x = -320; x < 620; x += 10) {
B_x = B_y = B_z = 0;
for (ptr = segments; ptr; ptr = ptr->next) {
a_1 = ptr->start.x;
a_2 = ptr->start.y;
b_1 = ptr->end.x;
b_2 = ptr->end.y;
a_3 = b_3 = 0;
R2_A = SQR(x-a_1)+SQR(y-a_2)+SQR(z-a_3);
R2_B = SQR(x-b_1)+SQR(y-b_2)+SQR(z-b_3);
L2 = SQR(a_1-b_1)+SQR(a_2-b_2)+SQR(a_3-b_3);
alpha2 = acos ((L2 + R2_A - R2_B) / (2 * sqrt (L2) * sqrt (R2_A)));
shortest_distance = sqrt(R2_A*(1-SQR(cos(M_PI-alpha2))));
alpha1 = M_PI - asin (shortest_distance / sqrt(R2_B));
if (R2_A > R2_B)
if ((R2_A - SQR(shortest_distance)) > L2) /* then we are outside */
alpha1 = asin (shortest_distance / sqrt(R2_B));
term1 = current * mu_naught / (4*M_PI*shortest_distance);
term2 = cos (alpha2);
term3 = cos (alpha1);
the_gnu_who_knew = term1 * (term2 - term3);
x_hat = (b_2 - a_2) * (z - a_3) - (b_3 - a_3) * (y - a_3);
y_hat = (b_3 - a_3) * (x - a_1) - (b_1 - a_1) * (z - a_3);
z_hat = (b_1 - a_1) * (y - a_2) - (b_2 - a_2) * (x - a_1);
cat_in_the_hat = sqrt(SQR(x_hat)+SQR(y_hat)+SQR(z_hat));
x_hat /= cat_in_the_hat;
y_hat /= cat_in_the_hat;
z_hat /= cat_in_the_hat;
B_x += x_hat * the_gnu_who_knew;
B_y += y_hat * the_gnu_who_knew;
B_z += z_hat * the_gnu_who_knew;
}
printf ("%g\n", sqrt (SQR(B_x)+SQR(B_y)+SQR(B_z)));
}
printf
("\n");
}
}
int main (int argc, char
**argv)
{
char buf[256];
int begin = 1;
struct point
first, current, previous;
struct segment
*ptr;
double z;
if (argc != 2)
die
("Usage: %s z_height\n", argv[0]);
sscanf (argv[1],
"%lf", &z);
while (fgets
(buf, 256, stdin) != NULL) {
if
(sscanf (buf, "%lf %lf", ¤t.x, ¤t.y) == 2) {
if (begin) {
first = current;
begin = 0;
} else {
add_segment (previous, current);
}
previous = current;
}
else {
fprintf (stderr, "empty line ... starting a new loop\n");
add_segment (previous, first);
begin = 1;
}
}
for (ptr = segments;
ptr; ptr = ptr->next)
fprintf
(stderr, "%g\t%g\n%g\t%g\n\n", ptr->start.x, ptr->start.y, ptr->end.x,
ptr->end.y);
field (z);
return 0;
}
Programs for the Experimental
Procedure
The following program was written to interface the xy translational device with the Labmaster. It moves the SMP in steps along the x and y directions. This program also triggers the TMS coil and collects the data for the measurements. The command line to run this program is:
control lowX hiX lowY hiY stepX stepY
The following is the actual code for the program.
#include <labmaster.h>
#define LPT 0x378
#define out_dac(a,b)
{int z; for (z = 0; z < 999; z++); out_dac (a, b);}
int main (int argc,
char **argv)
{
int x, y;
int lowX, hiX,
lowY, hiY, stepX, stepY;
int dummy1,dummy2,dummy3,
reading1,reading2, reading3;
FILE *gnuplot1,
*gnuplot2, *gnuplot3, *triplet;
gnuplot1 = fopen
("blah1.dat", "w");
gnuplot2 = fopen
("blah2.dat", "w");
gnuplot3 = fopen
("blah3.dat", "w");
triplet = fopen
("blahtrp.dat", "w");
if ((gnuplot1
== NULL) || (triplet == NULL)) {
fprintf
(stderr, "Facist bastard one won't save to disk!\n");
exit
(1);
}
if ((gnuplot2
== NULL) || (triplet == NULL)) {
fprintf
(stderr, "Facist bastard two won't save to disk!\n");
exit
(1);
}
if ((gnuplot3
== NULL) || (triplet == NULL)) {
fprintf
(stderr, "Facist bastard three won't save to disk!\n");
exit
(1);
}
labmaster_initialize
();
if (argc != 7)
{
fprintf
(stderr, "You suck!\nprogname lowX hiX lowY hiY stepX stepY\n");
exit
(1);
}
sscanf (argv[1],
"%d", &lowX);
sscanf (argv[2],
"%d", &hiX);
sscanf (argv[3],
"%d", &lowY);
sscanf (argv[4],
"%d", &hiY);
sscanf (argv[5],
"%d", &stepX);
sscanf (argv[6],
"%d", &stepY);
set_port_mode
(PORT_A,0);
printf ("Scanning
x from %d to %d in steps of %d\n",
lowX, hiX, stepX);
printf ("Scanning
y from %d to %d in steps of %d\n",
lowY, hiY, stepY);
for (x = lowX;
x < hiX; x += stepX) { //control x plane
out_dac
(0, x);
for
(y = lowY; y < hiY; y+=stepY) { //control y plane
out_dac (1, y);
delay(1000);
in_adc (1);
in_adc (1);
in_adc (2);
in_adc (2);
in_adc (3);
in_adc (3);
write_port(PORT_A,-1);
delay (60);
in_adc (1);
dummy1 = in_adc (1);
in_adc (2);
dummy2 = in_adc (2);
in_adc (3);
dummy3 = in_adc (3);
write_port(PORT_A,0);
delay (60);
write_port(PORT_A,-1);
outp (LPT, 0xff);
delay (60);
outp (LPT, 0);
in_adc (1);
reading1 = in_adc (1);
in_adc (2);
reading2 = in_adc (2);
in_adc (3);
reading3 = in_adc (3);
write_port(PORT_A,0);
printf ("%d %d %d %d %d\n", x, y,dummy1-reading1, dummy2-reading2,dummy3-reading3);
fprintf (triplet, "%d %d %d %d %d\n", x, y,dummy1-reading1,dummy2-reading2,
dummy3-reading3);
fprintf (gnuplot1, "%d \n",dummy1-reading1);
fprintf (gnuplot2, "%d \n",dummy2-reading2);
fprintf (gnuplot3, "%d \n",dummy3-reading3);
}
fprintf (gnuplot1, "\n");
fprintf (gnuplot2, "\n");
fprintf (gnuplot3, "\n");
}
fclose (triplet);
fclose (gnuplot1);
fclose (gnuplot2);
fclose (gnuplot3);
return 0;
}