path: root/fiz/naloga/numerično.c

#include <stdlib.h>
#include <stdio.h>
#include <error.h>
#include <math.h>
#include <signal.h>
#include <string.h>
#include <errno.h>
#define UVOD "program za numerični izračun jakosti magnetnega polja okoli helmholtzove tuljave\n" \
	"sem spisal anton luka šijanec za projektno nalogo pri fiziki v tretjem letniku gimb.\n" \
	"uporaba: %s	in nujni argumenti po vrsti:\n" \
	"	1. radij enega navitja v metrih\n" \
	"	2. tok, ki teče po vodniku v amperih\n" \
	"	3. število navojev na enem navitju\n" \
	"	4. razmak med merilnimi točkami v metrih\n" \
	"	5. koliko meritev od središča v obe dimenziji naj napravimo\n" \
	"	6. koliko kotov naj ima navitje - računamo, kot da je mnogokotnik\n" \
	"nenujni argumenti: prazen niz argument nastavi na privzeto vrednost.\n" \
	"	7. tip izhodnih podatkov (pgm, ppm ali tsv). privzeto je to pgm.\n" \
	"	8. razmak med navitjima. privzeto je to r/2.\n" \
	"	9. zamik enega navitja v metrih. privzeto sta osi navitij ista premica.\n" \
	"	A. polmer druge tuljave. privzeto je enak polmeru prve tuljave.\n" \
	"	B. tok druge tuljave v metrih. privzeto je enak toku prve tuljave.\n" \
	"oblika izhodnih podatkov, če je 7. parameter pgm, so pgm slike z vrednostmi 0-255\n" \
	"	- slika je prerez tuljave. magnetno polje teče vodoravno.\n" \
	"	- vrednosti direktno korelirajo z izračunano jakostjo v decigaussih: 10e-5 tesla\n" \
	"	- slika je široka 1+2*koliko in visoka 1+2*koliko (5. argument) slikovnih točk\n" \
	"oblika izhodnih podatkov, če je 7. parameter ppm, so barvne slike z 8 biti na barvo\n" \
	"	- barvi rdeča in zelena predstavljata komponenti vektorjev polja i in j\n" \
	"oblika izhodnih podatkov, če je 7. parameter tsv, je tsv, z naslednjimi stolpci:\n" \
	"	1. komponenta i krajevnega vektorja točke meritve. 0,0 je v središču tuljave.\n" \
	"	2. komponenta j krajevnega vektorja točke meritve. 0,0 je v središču tuljave.\n" \
	"	3. komponenta i vektorja magnetnega polja v točki meritve.\n" \
	"	4. komponenta j vektorja magnetnega polja v točki meritve.\n" \
	"	5. jakost magnetnega polja v teslah - tokrat ni v decigaussih!\n" \
	"način izdelave animacije: podan naj bo samo en argument - animacija - TOLE NE DELA\n" \
	"	- delajo se datoteke animacija0000.ppm, animacija0001.ppm, ...\n" \
	"	- parametri se sinusoidno spreminjajo po vdelanih konstantah.\n" \
	"	- slike lahko recimo s ffmpeg(1) nato pretvorite v videoposnetek."
enum oblika {
	PGM,
	PPM,
	TSV,
	ANIMACIJA
};
struct vektor {
	long double i; // x - desno na sliki
	long double j; // y - gor na sliki
	long double k; // z - v monitor
};
struct vektor seštej (struct vektor a, struct vektor b) {
	struct vektor r = {
		.i = a.i + b.i,
		.j = a.j + b.j,
		.k = a.k + b.k,
	};
	return r;
}
struct vektor vektorski_produkt (struct vektor a, struct vektor b) {	// ne bom implementiral
	struct vektor r = {						// matrik
		.i = a.j*b.k - a.k*b.j,
		.j = a.k*b.i - a.i*b.k,
		.k = a.i*b.j - a.j*b.i
	};
	return r;
}
struct vektor množi (struct vektor a, long double d) {
	struct vektor r = {
		.i = a.i * d,
		.j = a.j * d,
		.k = a.k * d
	};
	return r;
}
long double absolutno (struct vektor a) {
	return sqrtl(a.i*a.i+a.j*a.j+a.k*a.k);
}
#define MU0 4e-6*M_PI
struct vektor tuljava (long double R,	unsigned kotov, struct vektor m /* meritev - krajevni */) {
	long double dl_abs = 2*M_PI*R/kotov;	// metri - dolžina vodnika
	long double dr = 2*M_PI/kotov;		// radiani - kot med dl in točko na (0,R - vrh zanke)
	struct vektor B = {
		.i = 0,
		.j = 0,
		.k = 0
	};
	for (unsigned i = 0; i < kotov; i++) {
		long double theta = dr*i;	// kot na krogu
		struct vektor dl;
		dl.j = cosl(theta)*dl_abs;
		dl.k = -sinl(theta)*dl_abs;	// minus po skici sodeč ://of.sijanec.eu/sfu/skic.jpg
		dl.i = 0;			// sicer je vseeno, m je na z = 0 in gledamo vse
		struct vektor r;
		r.i = 0;
		r.j = sinl(theta)*R;
		r.k = cosl(theta)*R;
		r = seštej(r, m);
		B =	seštej(B,
				množi(
					vektorski_produkt(dl, r),
					1/(absolutno(r)*absolutno(r)*absolutno(r))
				)
			);
	}
	B = množi(B, MU0/(4*M_PI));
	return B;
} // ena zanka ob toku 1 A. pomnoži s tokom in številom navitij. 0,0 je v sredini. B teče v desno.
struct vektor zavrti_okoli (struct vektor točka, struct vektor središče, long double kot) {
	točka = seštej(točka, množi(središče, -1));
	struct vektor r = {
		.i = cosl(kot)*točka.i + sin(kot)*točka.j,
		.j = cosl(kot)*točka.j - sin(kot)*točka.i
	};
	return seštej(r, središče);
} // TODO: implement
void natisni (FILE * f, struct vektor v, const char * i) {
	fprintf(f, "vektor %s {\n\t.i = %Lf,\n\t.j = %Lf,\n\t.k = %Lf\n}\n", i, v.i, v.j, v.k);
}
int nariši (long double R, long double I, unsigned n, long double razmak, int koliko,
		unsigned kotov, enum oblika oblika, long double med_tuljavama,
		struct vektor zamik_izven_osi, long double R2, long double I2, FILE * out) {
	struct vektor merilno_mesto = {	// krajevni vektor
		.k = 0
	};
	if (oblika == PGM)
		fprintf(out, "P5 %u %u 255\n", koliko*2+1, koliko*2+1);
	if (oblika == PPM)
		fprintf(out, "P6 %u %u 255\n", koliko*2+1, koliko*2+1);
	struct vektor Rpolovic = {
		.i = med_tuljavama,
		.j = 0,
		.k = 0
	};
	for (int i = -koliko; i <= koliko; i++) {
		merilno_mesto.i = i*razmak;
		for (int j = -koliko; j <= koliko; j++) {
			merilno_mesto.j = j*razmak;
			struct vektor merilno_mesto2 = seštej(
					merilno_mesto,
					zamik_izven_osi
				);
			struct vektor B =
				seštej(
					množi(
						množi(
							tuljava(
								R,
								kotov,
								seštej(
									merilno_mesto,
									Rpolovic
								)
							),
							n
						),
						I
					),
					množi(
						množi(
							tuljava(
								R2,
								kotov,
								seštej(
									merilno_mesto2,
									množi(
										Rpolovic,
										-1
									)
								)
							),
							n
						),
						I2
					)
				);
			switch (oblika) {
				case PGM:
					if (10000*absolutno(B) > 255)
						putc(255, out);
					else
						putc(10000*absolutno(B), out);
					break;
				case PPM:
#define NATISNI_KOMPONENTO(x)		if (10000*fabsl(B.x) > 255) \
						putc(255, out); \
					else \
						putc(10000*fabsl(B.x), out);
					NATISNI_KOMPONENTO(i);
					NATISNI_KOMPONENTO(j);
					NATISNI_KOMPONENTO(k);
					break;
				case TSV:
					fprintf(out, "%Lf\t%Lf\t%Lf\t%Lf\t%Lf\n",
							merilno_mesto.i, merilno_mesto.j,
							B.i, B.j, absolutno(B));
					break;
				case ANIMACIJA:
					abort(); // invalid
					break;
			}
		}
	}
	return 0;
}
int shouldexit = 0;
void handler (int s __attribute__((unused))) {
	shouldexit++;
}
int main (int argc, char ** argv) {
	if (argv[1] && !strcmp(argv[1], "animacija"))
		goto animacija;
	if (argc < 8)
		error(1, 0, UVOD, argv[0] ? argv[0] : "./numerično");
	long double R = strtold(argv[1], NULL);
	long double I = strtold(argv[2], NULL);
	unsigned n = strtol(argv[3], NULL, 10);
	long double razmak = strtold(argv[4], NULL);
	int koliko = strtold(argv[5], NULL);
	unsigned kotov = strtol(argv[6], NULL, 10);
	enum oblika oblika = PGM;
	long double med_tuljavama = R/2;
	struct vektor zamik_izven_osi = {0, 0, 0};
	long double R2 = R;
	long double I2 = I;
	if (argc > 7 && argv[7][0])
		switch (argv[7][1]) {
			case 'G':
			case 'g':
				oblika = PGM;
				break;
			case 'P':
			case 'p':
				oblika = PPM;
				break;
			case 'S':
			case 's':
				oblika = TSV;
				break;
			case 'N':
			case 'n':
				oblika = ANIMACIJA;
				break;
		}
	if (argc > 8 && argv[8][0])
		med_tuljavama = strtold(argv[8], NULL);
	if (argc > 9 && argv[9][0])
		zamik_izven_osi.j = strtold(argv[9], NULL);
	if (argc > 0xA && argv[0xA][0])
		R2 = strtold(argv[0xA], NULL);
	if (argc > 0xB && argv[0xB][0])
		I2 = strtold(argv[0xB], NULL);
	if (oblika != ANIMACIJA)
		return nariši(R, I, n, razmak, koliko, kotov, oblika, med_tuljavama,
			       	zamik_izven_osi, R2, I2, stdout);
animacija:
	signal(SIGTERM, handler);
	signal(SIGINT, handler);
	unsigned čas = 0;
	while (!shouldexit && čas < 21) {
		char fn[25] = "animacija";
		sprintf(fn+strlen(fn), "%04d.ppm", čas);
		fprintf(stderr, "%s\n", fn);
		FILE * f = fopen(fn, "w");
		if (!f)
			error(2, errno, "fopen");
		long double R = 0.06;
		long double I = 2;
		unsigned n = 320;
		long double razmak = 0.0007;
		unsigned koliko = 350;
		unsigned kotov = 30;
		enum oblika oblika = PPM;
		long double med_tuljavama = (R/2/10)*(čas);
		struct vektor zamik = {
			.i = 0,
			.j = 0, // zamik osi
			.k = 0
		};
		long double R2 = R;
		long double I2 = I;
		if (nariši(R, I, n, razmak, koliko, kotov, oblika, med_tuljavama, zamik, R2, I2, f))
			error(3, errno, "nariši");
		if (fclose(f))
			error(4, errno, "fclose");
		čas++;
	}
	return 0;
}