244 lines
6.5 KiB
C
244 lines
6.5 KiB
C
/* Copyright (C) 2013 Riccardo Greco rigreco.grc@gmail.com.
|
|
*
|
|
* This project is based on 1999-2000 Thesis work of Greco Riccardo.
|
|
* It implement an Runge Kutta 4(5)^ order integration numerical method of differential equations set
|
|
* by use of double precision floating point operation in Aztec C65 language.
|
|
* It allow to simulate different mathematical models such as:
|
|
* Resistance Capacitor electrical circuit, Direct Current electric motor,
|
|
* Alternative Current three phase induction motor.
|
|
*
|
|
* Thanks to Bill Buckels for his invaluable support:
|
|
* Aztec C compilers http://www.aztecmuseum.ca/compilers.htm
|
|
*/
|
|
/* PROCESS MODULE */
|
|
#include <math.h>
|
|
#include <stdio.h>
|
|
#include <fcntl.h>
|
|
#include <errno.h>
|
|
#define chk 1 /* set chk to 1 charge */
|
|
#define rkfr 1.0e-6 /* RKF max resolution h<rk4r */
|
|
/*#define hmax 0.25 /* need for adaptative step
|
|
#define emax 20.0e-3
|
|
#define emin 10.0e-3*/
|
|
|
|
/* GENERAL GLOBAL */
|
|
extern double sur,r,c,tau,h,time,time2,yinit,keep;
|
|
extern char fbuf[80];
|
|
extern int cnt,fd;
|
|
extern struct data {
|
|
char v[20];
|
|
char i[20];
|
|
char ideal[20];
|
|
};
|
|
|
|
/* LOCAL GLOBAL */
|
|
char *outr[4]={"t="," vc(t)="," vcerr(t)="," i(t)="};
|
|
char *FILE_NAME="rcrk";
|
|
double swi;
|
|
|
|
/*Fehlberg parameters */
|
|
double acca[5]={1.0/4.0,3.0/8.0,12.0/13.0,1.0,1.0/2.0};
|
|
double erre1[5]={1.0/4.0,0.0,0.0,0.0,0.0};
|
|
double erre2[5]={3.0/32.0,9.0/32.0,0.0,0.0,0.0};
|
|
double erre3[5]={1932.0/2197.0,-7200.0/2197.0,7296.0/2197.0,0.0,0.0};
|
|
double erre4[5]={439.0/216.0,-8.0,3680.0/513.0,-845.0/4104.0,0.0};
|
|
double erre5[5]={-8.0/27.0,2.0,-3544.0/2565.0,1859.0/4104.0,-11.0/40.0};
|
|
|
|
double erre[5]={1.0/360.0,-128.0/4275.0,-2197.0/75240.0,1.0/50.0,2.0/55.0};
|
|
|
|
double ips5[5]={25.0/216.0,0.0,1408.0/2565.0,2197.0/4104.0,-1.0/5.0};
|
|
/*double ips6[6]={16.0/135.0,0.0,6656.0/12825.0,28561.0/56430.0,-9.0/50.0,2.0/55.0};*/
|
|
|
|
|
|
|
|
/* Functions sets */
|
|
int funrk(x,y,frk)
|
|
double x,y,*frk;
|
|
{
|
|
frk[0] = (double)(((sur*swi)-y)/tau); /* RC Circuit Voltage frk[0]=vc'(t) Charge differential equation */
|
|
} /* vc'(t)=(1/RC)*(E-vc(t)) charge*/
|
|
/* vc'(t)=(1/RC)*(-vc(t)) discharge */
|
|
double current(v)
|
|
double(v);
|
|
{
|
|
return(((sur*swi)-v)/r); /* i(t)=(E-vc(t))/R charge*/
|
|
} /* i(t)=(-vc(t))/R discharge */
|
|
|
|
double exact(x)
|
|
double(x);
|
|
{
|
|
return (((yinit-sur)*exp(-x/tau))+sur); /* Exact solution charge vc(t)=(vc(0)-E)*exp(-t/RC)+E */
|
|
}
|
|
|
|
double exact2(x)
|
|
double (x);
|
|
{
|
|
return (keep*exp(-x/tau)); /* Exact solution discharge vc(t)=(vc(td))*exp(-t/RC) */
|
|
}
|
|
|
|
|
|
|
|
/*RKF45 Module */
|
|
int runge_kutta(x,y,h,y3)
|
|
double x,y,h,*y3;
|
|
{
|
|
double r1,r2,r3,r4,r5,r6,frk,err; /*sigma,div,div2;
|
|
|
|
/* RK inizialization */
|
|
y3[0] = (double)0.0;
|
|
frk = (double)0.0;
|
|
err=(double)0.0;
|
|
|
|
/* do
|
|
{
|
|
|
|
/* 1^ order */
|
|
funrk(x,y,&frk);
|
|
r1=frk;
|
|
/* 2^ order */
|
|
funrk(x+acca[0]*h,(double)y+erre1[0]*r1*h,&frk);
|
|
r2=frk;
|
|
/* 3^ order */
|
|
funrk(x+acca[1]*h,(double)y+erre2[0]*r1*h+erre2[1]*r2*h,&frk);
|
|
r3=frk;
|
|
/* 4^ order */
|
|
funrk(x+acca[2]*h,(double)y+erre3[0]*r1*h+erre3[1]*r2*h+erre3[2]*r3*h,&frk);
|
|
r4=frk;
|
|
/* 5^ order */
|
|
funrk(x+acca[3]*h,(double)y+erre4[0]*r1*h+erre4[1]*r2*h+erre4[2]*r3*h+erre4[3]*r4*h,&frk);
|
|
r5=h*frk;
|
|
|
|
/* 6^ order only for estimated error for adaptative step method
|
|
funrk(x+acca[4]*h,(double)y+erre5[0]*r1*h+erre5[1]*r2*h+erre5[2]*r3*h+erre5[3]*r4*h+erre5[4]*r5*h,&frk);
|
|
r6=frk;
|
|
|
|
err=(fabs(erre[0]*r1+erre[1]*r3+erre[2]*r4+erre[3]*r5+erre[4]*r6)*h); /* estimate the error
|
|
|
|
h=h/2; /*non adaptative*/
|
|
|
|
/* adaptative step integration
|
|
div=(double)(rk4r/err);
|
|
div2=(double)pow(div,0.25);
|
|
sigma=(double)(0.84*div2);
|
|
|
|
if (sigma<=0.1) h=0.1*h;
|
|
else if (sigma>=4) h=4*h;
|
|
else h=sigma*h;
|
|
/* Alternative method
|
|
if (err>=emax) h=h/2;
|
|
if (err<=emin) h=h*2;
|
|
|
|
if (h>hmax) h=hmax; */
|
|
/*
|
|
}
|
|
while (err>rkfr);
|
|
*/
|
|
|
|
/* y3 weighted average of operators 5^ order */
|
|
y3[0]=(double)(y+(ips5[0]*r1+ips5[2]*r3+ips5[3]*r4+ips5[4]*r5)*h);
|
|
}
|
|
|
|
ovmain()
|
|
{
|
|
char buf[20],rev;
|
|
double atof(); /* DON'T FORGET THIS */
|
|
double x,x2,nc,y3,y,errtx,cur,exa;
|
|
struct data rk;
|
|
|
|
|
|
/* inizialization */
|
|
cnt=(int)0.0; /* Reset counter */
|
|
y3 = (double)0.0;
|
|
x=(double)0.0; /* time in simulation */
|
|
x2=(double)0.0; /* time in discharge */
|
|
keep=(double)0.0;
|
|
|
|
swi=(double)1.0; /* switch set to 1 = charge */
|
|
y=yinit; /* set initial condition parameter vc(o) */
|
|
/*h=hmax; /* need for adaptative step */
|
|
|
|
/* Start simulation */
|
|
scr_clear();
|
|
scr_curs(0,0);
|
|
puts("Pocessing...");
|
|
|
|
/* open and write header data to file */
|
|
open(FILE_NAME,O_WRONLY|O_APPEND,0xC3);
|
|
write(fd,fbuf,80);
|
|
|
|
buf[0]='\0'; /* buffers reset */
|
|
fbuf[0]='\0';
|
|
|
|
/* RK main Cycle */
|
|
while (x<time) {
|
|
|
|
runge_kutta(x,y,h,&y3);
|
|
cur=current(y3);
|
|
|
|
if (x <= time2)
|
|
{
|
|
exa=(double)exact(x);
|
|
errtx=(double)fabs(y3-exa); /* exact solution in charge */
|
|
}
|
|
|
|
if (x > time2)
|
|
{
|
|
if (chk)
|
|
{
|
|
swi=(double)0.0; /* switch set to 0 = discharge */
|
|
keep=(double)y3; /* keep the value of vc(td) in the istant of discharge */
|
|
!chk; /* NOT chk */
|
|
}
|
|
exa=(double)exact2(x2);
|
|
errtx=(double)fabs(y3-exa); /* exact solution in discharge */
|
|
x2=x2+h; /* increment time in discharge */
|
|
}
|
|
|
|
|
|
/* store output in buffer */
|
|
ftoa(x,buf,6,2);
|
|
strcat(fbuf,outr[0]);
|
|
strcat(fbuf,buf); /* store time (x) */
|
|
|
|
ftoa(y3,buf,6,2);
|
|
strcat(fbuf,outr[1]);
|
|
strcat(fbuf,buf); /* store voltage vc(t) (y3) */
|
|
|
|
ftoa(errtx,buf,6,2);
|
|
strcat(fbuf,outr[2]);
|
|
strcat(fbuf,buf); /* store voltage error vcerr(t) (errtx) */
|
|
|
|
ftoa(cur,buf,6,2);
|
|
strcat(fbuf,outr[3]);
|
|
strcat(fbuf,buf); /* store current i(t) (cur) */
|
|
|
|
/* converting using max floating resolution */
|
|
ftoa(y3,rk.v,14,1); /* vc(t) */
|
|
ftoa(cur,rk.i,14,1); /* i(t) */
|
|
ftoa(exa,rk.ideal,14,1); /* vc(t) exact solution */
|
|
|
|
puts(fbuf); /* output to screen */
|
|
|
|
/* write data result to file */
|
|
write(fd,rk.v,20); /* write vc(t) */
|
|
write(fd,rk.i,20); /* write i(t) */
|
|
write(fd,rk.ideal,20); /* write vc(t) exact solution */
|
|
|
|
fbuf[0]='\0'; /* reset buffer */
|
|
|
|
y=y3;
|
|
x=x+h; /* increment time in simulation */
|
|
cnt++; /* renew counter */
|
|
|
|
}
|
|
close(fd); /* close the data file */
|
|
|
|
/*printf("Cycles number= %d",(cnt-1)); /* NOT YET itoa ?*/
|
|
|
|
puts("Press any key to start plotting...");
|
|
getch();
|
|
|
|
return 0;
|
|
}
|
|
|