ModSim: Almost finished implementing models.

modsim/ass3/Makefile

@@ -8,6 +8,6 @@ PROGS=test harm_osc
 all: ${PROGS}
 
 test: test.o logger.o methods.o
-harm_osc: harm_osc.o logger.o methods.o
+harm_osc: harm_osc.o models.o logger.o methods.o
 
 clean: ; rm -vf *.o ${PROGS}

modsim/ass3/harm_osc.c

@@ -2,80 +2,50 @@
 #include <stdio.h>
 #include <math.h>
 #include "methods.h"
+#include "models.h"
 #include "logger.h"
 
-#define DT .001
+int main(int argc, char **argv) {
+    if( argc < 2 ) {
+        printf("Usage: %s ODE\n", argv[0]);
+        return EXIT_FAILURE;
+    }
 
-#define INTEGRATE(method, t0, t1, FN) { \
-    method(t0, t1, DT, y0, y1, 2, &f_ ## FN, args_ ## FN); \
-    /*if( method(t0, t1, DT, y0, y1, 2, &f_ ## FN, args_ ## FN) ) puts("FAIL");*/ \
-    /* else printf("y1 using %-12s s:%.10f\n", #method ":", *y1); */\
-}
-
-#define COMPARE(t0, t1, fn) { \
-    /*printf("Function %s\n", #fn);*/ \
-    INTEGRATE(Euler, t0, t1, fn); \
-    INTEGRATE(RungeKutta2, t0, t1, fn); \
-    INTEGRATE(RungeKutta4, t0, t1, fn); \
-}
-
-#define PARAM(n) ((double *) params)[n]
-
-int f_osc(double t, double *y, double *dy, void *params) {
-    // s' = v
-    dy[0] = y[1];
-    // v' = -k * s
-    dy[1] = - PARAM(0) * y[0];
-    return 0;
-}
-
-int f_forced_osc(double t, double *y, double *dy, void *params) {
-    // s' = v
-    dy[0] = y[1];
-    // v' = -k * s - D * v + sin(w * t)
-    dy[1] = - PARAM(0) * y[0] - PARAM(1) * y[1] + sin(PARAM(2) * t);
-    return 0;
-}
-
-int f_lottka_volt(double t, double *y, double *dy, void *params) {
-    // x' = -a * x + c * d * x * y
-    dy[0] = - PARAM(0) * y[0] + PARAM(2) * PARAM(3) * y[0] * y[1];
-    // y' = b * y - d * x * y
-    dy[1] = PARAM(1) * y[1] - PARAM(3) * y[0] * y[1];
-    return 0;
-}
-
-int f_modified_pp(double t, double *y, double *dy, void *params) {
-    // x' = -a * x + c * d * x * y
-    dy[0] = - PARAM(0) * y[0] + PARAM(2) * PARAM(3) * y[0] * y[1];
-    // y' = b * ( 1 - y / y_m ) * y - d * x * y
-    //    = b * y - b * y * y / y_m - d * x * y
-    double by = PARAM(1) * y[1];
-    dy[1] = by - by * y[1] / PARAM(4) - PARAM(3) * y[0] * y[1];
-    return 0;
-}
-
-int main(void) {
     double y0[2] = {.0, 5.0},
            y1[2] = {.0, .0};
 
-    double args_osc[1]         = {.5},                     // {k}
+    const double args_osc[1]   = {.5},                     // {k}
            args_forced_osc[3]  = {1.0, 0.1, 2.0},          // {k, D, w}
            args_lottka_volt[4] = {.5, 1.0, .1, .1},        // {a, b, c, d}
-           args_modified_pp[5] = {.5, 1.0, .1, .1, 100.0}; // {a, b, c, d, y_m}
+           args_modified_pp[5] = {.5, 1.0, .1, .1, 100.0}, // {a, b, c, d, y_m}
+           args_pid_control[5] = {-2.0, -1.0, .1, 10.0, .001}; // {P, D, a, M, d}
     //printf("number of args: %d\n", sizeof(args_osc) / sizeof(double));
-    COMPARE(.0, 10.0, osc);
-    COMPARE(.0, 10.0, forced_osc);
-
-    // TODO: this model is not yet completed.
-    //printf("Function: lottka_volt\n");
-    INTEGRATE(RungeKutta4, .0, 10.0, lottka_volt);
 
-    // TODO: this model is not yet completed.
-    //printf("Function: modified_pp\n");
-    INTEGRATE(RungeKutta4, .0, 10.0, modified_pp);
+    switch( atoi(argv[1]) ) {
+        case 0: COMPARE(.0, 10.0, osc); break;
+        case 1: INTEGRATE(RungeKutta4, .0, 20.0, forced_osc); break;
+        case 2: INTEGRATE(RungeKutta4, .0, 20.0, lottka_volt); break;
+        case 3: INTEGRATE(RungeKutta4, .0, 20.0, modified_pp); break;
+        case 4:
+            y0[0] = .0;  // = v(0)
+            y0[1] = 3.0; // = s(0)
+            buf_size = ceil(args_pid_control[4] / DT);
+            v_buf = malloc(buf_size * sizeof(double));
+            s_buf = malloc(buf_size * sizeof(double));
+
+            for( int i = 0; i < buf_size; i++ ) {
+                v_buf[i] = y0[0];
+                s_buf[i] = y0[1];
+            }
+
+            INTEGRATE(RungeKutta4, .0, 100.0, pid_control);
+            free(v_buf);
+            free(s_buf);
+        break;
+        default: puts("ODE does not exist!"); return EXIT_FAILURE;
+    }
 
     logger_close();
 
-    return 0;
+    return EXIT_SUCCESS;
 }

modsim/ass3/logger.c

@@ -13,18 +13,13 @@ void logger(double t, int N, double *y) {
         logfile = stdout;
 
     fprintf(logfile, "%e", t);
-    if( N > 0 ) {
-        for( int i = 0; i < N - 1; i++ )
-            fprintf(logfile, ",%e", y[0]);
-        fprintf(logfile, ",%e\n", y[N-1]);
-    }
-    else
+    for( int i = 0; i < N; i++ )
+        fprintf(logfile, ",%e", y[i]);
     fprintf(logfile, "\n");
 }
 
 void logger_close() {
-    if( logfile && fclose(logfile) ) {
+    if( logfile && logfile != stdout && fclose(logfile) )
         printf("Failed to close logfile. Error no: %d", errno);
-    }
 }
 

modsim/ass3/models.c

+#include <math.h>
+#include "models.h"
+
+int pos = 0;
+
+int f_osc(double t, double *y, double *dy, void *params) {
+    // s' = v
+    dy[0] = y[1];
+    // v' = -k * s
+    dy[1] = - PARAM(0) * y[0];
+    return 0;
+}
+
+int f_forced_osc(double t, double *y, double *dy, void *params) {
+    // s' = v
+    dy[0] = y[1];
+    // v' = -k * s - D * v + sin(w * t)
+    dy[1] = - PARAM(0) * y[0] - PARAM(1) * y[1] + sin(PARAM(2) * t);
+    return 0;
+}
+
+int f_lottka_volt(double t, double *y, double *dy, void *params) {
+    // x' = -a * x + c * d * x * y
+    dy[0] = - PARAM(0) * y[0] + PARAM(2) * PARAM(3) * y[0] * y[1];
+    // y' = b * y - d * x * y
+    dy[1] = PARAM(1) * y[1] - PARAM(3) * y[0] * y[1];
+    return 0;
+}
+
+int f_modified_pp(double t, double *y, double *dy, void *params) {
+    // x' = -a * x + c * d * x * y
+    dy[0] = - PARAM(0) * y[0] + PARAM(2) * PARAM(3) * y[0] * y[1];
+    // y' = b * ( 1 - y / y_m ) * y - d * x * y
+    //    = b * y - b * y * y / y_m - d * x * y
+    double by = PARAM(1) * y[1];
+    dy[1] = by - by * y[1] / PARAM(4) - PARAM(3) * y[0] * y[1];
+    return 0;
+}
+
+int f_pid_control(double t, double *y, double *dy, void *params) {
+    int buf_begin = (pos + 1) % buf_size;
+
+    // U(t) = P * s(t - d) + D * v(t - d)
+    double U = PARAM(0) * s_buf[buf_begin] + PARAM(1) * v_buf[buf_begin];
+    // F(t) = U(t) - a * v(t)
+    double F = U - PARAM(2) * y[0];
+    // v'(t) = F(t) / M
+    dy[0] = F / PARAM(3);
+    // s'(t) = v(t)
+    dy[1] = y[0];
+
+    v_buf[pos] = y[0];
+    s_buf[pos] = y[1];
+    pos = buf_begin;
+
+    return 0;
+}

modsim/ass3/models.h

+#ifndef __MODELS_H__
+#define __MODELS_H__
+
+#define DT .001
+
+#define INTEGRATE(method, t0, t1, FN) { \
+    method(t0, t1, DT, y0, y1, 2, &f_ ## FN, (void *) args_ ## FN); \
+    /*if( method(t0, t1, DT, y0, y1, 2, &f_ ## FN, args_ ## FN) ) puts("FAIL");*/ \
+    /* else printf("y1 using %-12s s:%.10f\n", #method ":", *y1); */\
+}
+
+#define COMPARE(t0, t1, fn) { \
+    /*printf("Function %s\n", #fn);*/ \
+    INTEGRATE(Euler, t0, t1, fn); \
+    INTEGRATE(RungeKutta2, t0, t1, fn); \
+    INTEGRATE(RungeKutta4, t0, t1, fn); \
+}
+
+#define PARAM(n) ((const double *) params)[n]
+
+double *v_buf, *s_buf;
+int buf_size;
+
+int f_osc(double t, double *y, double *dy, void *params);
+int f_forced_osc(double t, double *y, double *dy, void *params);
+int f_lottka_volt(double t, double *y, double *dy, void *params);
+int f_modified_pp(double t, double *y, double *dy, void *params);
+int f_pid_control(double t, double *y, double *dy, void *params);
+
+#endif // __MODELS_H__

modsim/ass3/plot.py

 #!/usr/bin/env python
-from sys import argv, exit, stdin
+from sys import argv, stdin
 from pylab import figure, clf, plot, show, savefig
 
 try:
@@ -21,7 +21,7 @@ for i in range(len(lines)):
 figure(1)
 clf()
 for i in range(len(y)):
-    plot(x, y[i], 'x')
+    plot(x, y[i], '-')
 if len(argv) == 3:
     savefig(argv[2])
 show()

telematica/ass1/async.py

 import asyncore, socket
-from Queue import Queue
 import re
 import random
 from asyncbase import AsyncBase, MAJOR_VERSION, MINOR_VERSION

telematica/ass1/doc/conf.py

@@ -174,7 +174,7 @@ htmlhelp_basename = 'Chatz0rdoc'
 # (source start file, target name, title, author, documentclass [howto/manual]).
 latex_documents = [
   ('index', 'Chatz0r.tex', u'Chatz0r Documentation',
-   u'Taddeüs Kroes, Sander van Veen', 'manual'),
+   u'Taddeüs Kroes (6054129), Sander van Veen (6167969)', 'manual'),
 ]
 
 # The name of an image file (relative to this directory) to place at the top of

telematica/ass1/server.py

@@ -3,7 +3,6 @@
 import asyncore
 import logging
 import logging.config
-import os
 import re
 import socket
 import sys
@@ -53,7 +52,7 @@ class Server(asyncore.dispatcher):
         try:
             conn, addr = self.accept()
             self.log.info('accepted client %s:%d' % addr)
-            client = self.connect_client(conn, addr)
+            self.connect_client(conn, addr)
         except socket.error:
             self.log.warning('warning: server accept() threw an exception.')
             return