The ultimate authority which determines the action of the REC-R operators and predicates is the C program rec-sero.c, whose verbatim listing follows:
#import "rec.h"
#import "sero.h"
#import "rectbl.h"
#import <math.h>
#import <dpsclient/psops.h>
#import <appkit/graphics.h>
# define HZAR 500 /* array length for graph horizons */
int options[4]={FALSE,FALSE,FALSE,FALSE};
extern double r_dblpar; /* value of REC's $xxx.yy$ */
int stackptr=-1;
int (*RECPot)()=&potl5;
int RECPotl[20] = {
&potl1, /* free particle */
&potl1, /* free particle */
&potl2, /* finite well */
&potl3, /* barrier */
&potl4, /* step */
&potl5, /* harmonic oscillator */
&potl6, /* quartic oscillator */
&potl7, /* factored dirac h.o. */
&potl8, /* gaussian dirac eqn */
&potl9, /* Mathieu equation */
&potl10, /* linear potential */
&potl11, /* coulomb potential */
&potl12, /* inverse square */
&potl13, /* damped harmonic osc */
&potl14, /* dirac harmonic osc */
&potl15, /* factored inv. square */
&potl15, /* factored inv. square */
&potl15, /* factored inv. square */
&potl18, /* dirac harmonic osc */
&potl19 /* dirac periodic potl */
};
int RECOffset = 0;
double yoffset = 0.0;
double RECEnergy = 0.5;
double RECMass = 1.0;
double RECWeight = 1.0;
double RECGPoint = 0.0;
double RECGStep = 0.1;
double RECKPoint = 0.0;
double RECKStep = 0.1;
double RECPtntl[4];
double RECSolun[4];
double RECPush[10];
double horhi[HZAR];
double horlo[HZAR];
double lastex, lastwy;
int horix=0;
extern int ptj;
extern double hi, em, egy;
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* definition of REC operators and predicates */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
int rprze() {return(options[0]==TRUE);}
int rpron() {return(options[1]==TRUE);}
int rprtw() {return(options[2]==TRUE);}
int rprth() {return(options[3]==TRUE);}
void ropla() {RECOffset+=1; horix=0;} /* increase offset, reset horizon */
int rprlb() {return RECPush[stackptr]>=0.0;} /* top of stack is positive */
void roplc() {RECPush[stackptr]=2.0*tanh(0.5*RECPush[stackptr]);}
void ropld() {RECKStep*=r_dblpar;} /* modify R-K step size */
void rople() {RECEnergy+=r_dblpar;} /* modify Energy */
void roplf() {} /* shorten step */
void roplg() {double ex, wy; /* graph additional points */
if (stackptr >= 10 || stackptr < 1) return;
wy=RECPush[stackptr--]+((double)RECOffset)*RECGStep+yoffset;
ex=RECPush[stackptr--]+((double)RECOffset)*RECGStep;
videoplot(ex,wy,1);
}
void roplh() {double x, y, z, ex, wy; int i; /* graph additional points */
if (horix==0) return;
if (stackptr >= 10 || stackptr < 1) return;
wy=RECPush[stackptr--]+((double)RECOffset)*RECGStep+yoffset;
ex=RECPush[stackptr--]+((double)RECOffset)*RECGStep;
i=horix+RECOffset;
if (wy>=horhi[i]) {
horhi[i]=wy;
if (lastwy>=horhi[i-1]) {
videoplot(lastex,lastwy,0);
videoplot(ex,wy,1);
} else {
z=(horhi[i]-horhi[i-1]-wy+lastwy);
if (z==0.0) x=lastex; else
x=lastex+(RECGStep*(lastwy-horhi[i-1]))/z;
y=lastwy+(x-lastex)*((wy-lastwy)/RECGStep);
videoplot(x,y,0);
videoplot(ex,wy,1);
}
}
if (wy<horhi[i] && lastwy>=horhi[i-1]) {
z=(horhi[i]-horhi[i-1]-wy+lastwy);
if (z==0.0) x=lastex; else
x=lastex+(RECGStep*(lastwy-horhi[i-1]))/z;
y=lastwy+(x-lastex)*((wy-lastwy)/RECGStep);
videoplot(lastex,lastwy,0);
videoplot(x,y,1);
}
if (wy<horlo[i]) {
horlo[i]=wy;
if (lastwy<horlo[i-1]) {
videoplot(lastex,lastwy,0);
videoplot(ex,wy,1);
} else {
z=(horlo[i]-horlo[i-1]-wy+lastwy);
if (z==0.0) x=lastex; else
x=lastex+(RECGStep*(lastwy-horlo[i-1]))/z;
y=lastwy+(x-lastex)*((wy-lastwy)/RECGStep);
videoplot(x,y,0);
videoplot(ex,wy,1);
}
}
if (wy>horlo[i] && lastwy<=horlo[i-1]) {
z=(horlo[i]-horlo[i-1]-wy+lastwy);
if (z==0.0) x=lastex; else
x=lastex+(RECGStep*(lastwy-horlo[i-1]))/z;
y=lastwy+(x-lastex)*((wy-lastwy)/RECGStep);
videoplot(lastex,lastwy,0);
videoplot(x,y,1);
}
lastex=ex; lastwy=wy;
if (horix<HZAR-1) horix++;
}
void ropli() {
stackptr++;
if (stackptr>=10) return;
RECPush[stackptr]=RECSolun[0];
} /* j[0][0] solution element */
void roplj() {
stackptr++;
if (stackptr>=10) return;
RECPush[stackptr]=RECSolun[1];
} /* [0][1] solution element */
void roplk() {
stackptr++;
if (stackptr>=10) return;
RECPush[stackptr]=RECSolun[2];
} /* [1][0] solution element */
void ropll() {
stackptr++;
if (stackptr>=10) return;
RECPush[stackptr]=RECSolun[3];
} /* [1][1] solution element */
void roplm() {RECMass+=r_dblpar; em=RECMass;}
void ropln() {
stackptr++;
if (stackptr>=10) return;
RECPush[stackptr]=r_dblpar;
} /* push the $$ register */
void roplo() {RECKPoint=r_dblpar;} /* set origin */
void roplp() {
stackptr++;
if (stackptr>=10) return;
RECPush[stackptr]=RECPush[stackptr-1];
} /* copy top element on stack */
void roplq() {}
void roplr() {serss(
RECSolun,
RECPtntl,
&RECKPoint,
RECKStep,
RECEnergy,
RECPot,
serr4);}
/* one step of fourth order Runge-Kutta integration */
void ropls() { /* push, then increment the step variable */
stackptr++;
if (stackptr>=10) return;
RECPush[stackptr]=RECGPoint;
RECGPoint+=RECGStep;
}
void roplt() {
stackptr++;
if (stackptr>=10) return;
RECPush[stackptr]=0.5*(RECSolun[0]+RECSolun[3]);
} /* half trace of solution matrix = cos phi */
void roplu() {serum(RECSolun);} /* put unit matrix as solution matrix */
void roplv() { /* push independent variable */
stackptr++;
if (stackptr>=10) return;
RECPush[stackptr]=RECKPoint;
}
void roplw() {RECWeight+=r_dblpar; hi=RECWeight;}
/* ---- upper case ----- */
void ropua() {int i; /* reset offset for plotted point */
horix=0;
RECOffset=0;
for (i=0; i<HZAR; i++) horhi[i]=-10.0;
for (i=0; i<HZAR; i++) horlo[i]= 10.0;
}
int rprub() {return (RECPush[stackptr] <= 1.0 && RECPush[stackptr] >= -1.0);}
/* top of stack lies between -1.0 and 1.0 */
void ropuc() {int c; /* set color */
PSstroke();
c = *r_pc++;
switch (c) {
case 'r': PSsetrgbcolor(1.0,0.0,0.0); break;
case 'g': PSsetrgbcolor(0.0,1.0,0.0); break;
case 'b': PSsetrgbcolor(0.0,0.0,1.0); break;
case 'o': PSsetrgbcolor(0.8,0.4,0.1); break;
case 'v': PSsetrgbcolor(1.0,0.0,1.0); break;
case 'c': PSsetcmykcolor(1.0,0.0,0.0,0.0); break;
case 'm': PSsetcmykcolor(0.0,1.0,0.0,0.0); break;
case 'y': PSsetcmykcolor(0.0,0.0,1.0,0.0); break;
case 'k': PSsetgray(NX_BLACK); break;
case 'd': PSsetgray(NX_DKGRAY); break;
case 'l': PSsetgray(NX_LTGRAY); break;
case 'w': PSsetgray(NX_WHITE); break;
default: PSsetgray(NX_BLACK); break;}
}
void ropud() {RECKStep=r_dblpar;} /* set R-K step size */
void ropue() {RECEnergy=r_dblpar;} /* set Energy */
int rpruf() {return (RECPtntl[1]*RECPtntl[2]>=0.0);} /* coeff parity */
void ropug() {double ex, wy; /* graph additional points */
if (stackptr >= 10 || stackptr < 1) return;
wy=RECPush[stackptr--]+((double)RECOffset)*RECGStep+yoffset;
ex=RECPush[stackptr--]+((double)RECOffset)*RECGStep;
videoplot(ex,wy,0);
}
void ropuh() {double ex, wy; int i; /* graph first point */
wy=RECPush[stackptr--]+((double)RECOffset)*RECGStep+yoffset;
ex=RECPush[stackptr--]+((double)RECOffset)*RECGStep;
lastex=ex; lastwy=wy;
i=horix+RECOffset;
if (wy>horhi[i]) horhi[i]=wy;
if (wy<horlo[i]) horlo[i]=wy;
if (horix<HZAR-1) horix++;
}
void ropui() {
stackptr++;
if (stackptr>=10) return;
RECPush[stackptr]=RECPtntl[0];
} /* j[0][0] solution element */
void ropuj() {
stackptr++;
if (stackptr>=10) return;
RECPush[stackptr]=RECPtntl[1];
} /* [0][1] solution element */
void ropuk() {
stackptr++;
if (stackptr>=10) return;
RECPush[stackptr]=RECPtntl[2];
} /* [1][0] solution element */
void ropul() {
stackptr++;
if (stackptr>=10) return;
RECPush[stackptr]=RECPtntl[3];
} /* [1][1] solution element */
void ropum() {RECMass=em=r_dblpar;}
void ropun() {}
void ropuo() { /* origin, increment - ind. var. */
if (stackptr >= 10 || stackptr < 1) return;
RECGStep=RECPush[stackptr--];
RECKPoint=RECPush[stackptr--];
}
void ropup() {if (stackptr<-1) {stackptr=-1; return;} else stackptr--;}
int rpruq() {return(RECPush[stackptr]*RECPush[stackptr]>25.0);}
/* test for minimum scale */
void ropur() {serr6(RECSolun,&RECKPoint,RECGStep,RECEnergy,RECPot);}
/* one step sixth order Runge-Kutta integration */
void ropus(){
if (stackptr >= 10 || stackptr < 1) return;
RECGStep=RECPush[stackptr--];
RECGPoint=RECPush[stackptr--];
} /* Initiallize the graphical step */
void roput() {
stackptr++;
if (stackptr>=10) return;
hi=RECWeight;
RECPush[stackptr]=serar(RECEnergy,RECPot,&serr6);
} /* half trace of solution matrix over a period */
void ropuu() {yoffset=r_dblpar;} /* offset y-origin */
void ropuv() {int i;
i=*r_pc++;
if (i<='9') i-='0'; else {if (i<='Z') i-='A'-10; else i-='a'-10;}
RECPot=RECPotl[i];
printf("RECPot = %2d\n",i);
} /* choose potential */
void ropuw() {RECWeight=hi=r_dblpar;} /* define weight */
/* ------------ */
void roppl() {RECKPoint+=RECKStep;} /* advance by one Runge-Kutta sized step */
void ropca() {}
void ropqm() { /* print reg values */
printf("z = %6.3f, E = %6.3f\n",RECKPoint,RECEnergy);
printf("Solution = (%6.3f,%6.3f,%6.3f,%6.3f)\n",
RECSolun[0],
RECSolun[1],
RECSolun[2],
RECSolun[3]);
printf("Potential = (%6.3f,%6.3f,%6.3f,%6.3f)\n",
RECPtntl[0],
RECPtntl[1],
RECPtntl[2],
RECPtntl[3]);
printf("\n");
}
int rprns() {return (int)(random()&01);} /* random 1/2 pred */
int rprps() {return random()%10==1;} /* random 1/10 pred */