// This is the c++ version of the Geometric Clutch Model version 2.
// It is called ModelC.cpp and was written by C.B.Lindemann in 1997.
// It was the version published in 1994 in Cell Motil. Cytoskeleton V29:141
// It was developed in 1992-1993 at Oakland University, Rochester MI.48309
// and is the intellectual property of Charles B. Lindemann.
// Developed with support from the National Science Foundation.

#include <iostream.h>
#include <graphics.h>
#include <fstream.h>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <ctype.h>
#include <math.h>
#include <time.h>

FILE *printer; // pointer for printer applications
// ****************************************************
// These are the declared global costants, variables and arrays
const float pi = 3.1415926;
int ROUTE=1;// intial and continuation routing
int PL, EX, N, Q, total=0;//plotting data from previous run
double sx[32];//position coordinates x and y
double sy[32];
double G[50] [31];//storage arrays for plotting a run
double K[50] [31];
double ang [31];// angles for 30 segments 1 to 30
double sang[31];//shear angles
double cur[31];// curvature values
double param[16];//MODELING PARAMETERS
char comm1[20];//globally accessible identifier for print outs
int graphdriver = DETECT, graphmode;// setup for graphics

// ****************************************************
// this programs specific prototypes are as follows:
void cls(void);
void load_angles(void);
void load_cilium(void);
void print_param(void);
double find_angle(double, double, double);
void run_model(void);
void graph_output(void);
// ****************************************************
main()
{
char menu;
char deflt;
start : ;
cls();
cout << "Welcome to the C version of the Geometric Clutch model \n";
cout << "Written by Dr. Charles B. Lindemann in 1997\n";
cout << "It is based on the 1993 version 7.1 in Qbasic.\n\n\n";
cout << "To run this program choose A it will automatically load the \n";
cout << "default settings for a 10micron cilium. \n";
cout << "IMPORTANT! CHOICES B and C will work only if you have a compatible \n";
cout << "printer and C requires that the printer takes PCL-5 OR PCL-6 GRAPHICS \n\n";
cout << "A) RUN MODEL.\n";
cout << "B) PRINT OUT MODELING PARAMETERS.\n";
cout << "C) PLOT LAST MODEL RUN.\n";
cout << "Q) EXIT FROM PROGRAM.\n\n\n";
cout << "\t Press the letter of the desired selection: ";
menu = toupper(getch());
cout << '\n';

switch (menu)
{
	case 'A':goto first;
	case 'B':goto second;
	case 'C':goto third;
	case 'Q': return 0;
	default :goto start;
}

first  : cls();
ROUTE = 1;
total = 0;
EX = 0;
load_cilium();
load_angles();
run_model();
goto start;
second  : cls();
print_param();
cout << "SENT TO PRINTER:"<< endl;
getch();
goto start;
third : cls();
graph_output();
goto start;
}
//clear screen function ************************
void cls(void)
{
// the next four lines actually clear the screen and return the cursor
int i=1;
for (i=1; i<26; i++)
{ cout << '\n';
}
gotoxy (1,1);// A function from conio.h that works dandy
return;
}
//**********************************************
// Load angles from storage and compute starting position
void load_angles(void)
{
 double Z = param[10];
 int i = 0;
for (i=1; i<31; i++) ang[i]=0;
// segment length Z is a modeling parameter
sx[1]=sx[0]=0.0;
sy[1]=sy[0]=0.0;
for (int j=1; j<31; j++) {
	sx[j+1] = sx[j] +(Z * cos(ang[j]));
	sy[j+1] = sy[j] +(Z * sin(ang[j]));
	cout << sx[j] <<"  " << sy[j] << "  " << ang[j] << endl;
	}
cur[0]=cur[1]=0.0;
sang[0]=sang[1]=0.0;
for ( i=2 ; i<31 ; i++){
	cur[i] = (ang[i]-ang[i-1]) / Z;//calculates the curvatures
	sang[i]=ang[i] - ang[1];
	cout << i << "  " << cur[i] << endl;
	}

}// end of load angles function
 //***************************************************
// Load default settings for a ten micron cilium
void load_cilium(void)
{
cout << "\nTHESE ARE THE DEFAULT SETTINGS FOR A 10 micron CILIUM: " << endl;
param[1]=.14;
param[2]=1.2e-07;
param[3]=.2;
param[4]=.06;
param[5]=.01;
param[6]=33000.;
param[7]=33000.;
param[8]=1.0e-13;
param[9]=1.0e-05;
param[10]=3.33e-05;
param[11]=0.028;
param[12]=7000.;
param[13]=.0001;
param[14]=1.0;
param[15]=1.0;
cout << "\tMODELING PARAMETERS" << endl;
cout << "\n";
cout << "A) COUPLING= "<< param[1] << endl;
cout << "B) BRIDGEFORCE= "<< param[2] << endl;
cout << "C) SNAP= "<< param[3] << endl;
cout << "D) SWITCH.P=" << param[4] << endl;
cout << "E) SWITCH.R=" << param[5] << endl;
cout << "F) HISTER.P=" << param[6] << endl;
cout << "G) HISTER.R=" << param[7] << endl;
cout << "H) IE= "<< param[8] << endl;
cout << "I) DIAM= " << param[9] << endl;
cout << "J) Z= " << param[10] << endl;
cout << "K) DRAG= " << param[11] << endl;
cout << "L) GAIN= " << param[12] << endl;
cout << "M) ITS= " << param[13] << endl;
cout << "N) STATE= " << param[14] << endl;
cout << "O) BASE= " << param[15] << endl;

}

//***********************************************
//Print model parameters function
void print_param(void)
{
// local variables
int i;
char menu;
char date[10];
char comm2[20];
cout << "ENTER THE DATE; ";
cin >> date;
cout << "OTHER INFORMATION:";
cin >> comm2;
printer =fopen ("prn", "w");
if (printer==0){
	 printf("printer was not located\n");
	 goto end;
	 }
fputs("\n\t ADVANCED MODEL version 7.1\n", printer);
fputs("\t   THE GEOMETRIC CLUTCH \n\n", printer);
fputs("   CGS UNITS and simplest torque balance\n\n", printer);
fprintf(printer,"  MODELING PARAMETERS     DATED: %s\n", date);
fprintf(printer,"  IDENTIFYING COMMENT: %s\n", comm1);
fprintf(printer,"  IDENTIFYING COMMENT: %s\n\n\n", comm2);
fprintf(printer,"SIDE TO SIDE COUPLING = %.2f",param[1]);
fprintf(printer,"\nFORCE PER DYNEIN BRIDGE IN DYNES(tug)= %.2e", param[2]);
fprintf(printer,"\nELASTIC CONST. OF PASSIVE ELASTIC NEXIN(snap)= %.2f", param[3]);
fprintf(printer,"\nSWITCH PROBABILITY P-BRIDGE = %.3f", param[4]);
fprintf(printer,"\nSWITCH PROBABILITY R-BRIDGE = %.3f", param[5]);
fprintf(printer,"\nCOEFFICIENT OF HISTERESIS(P) = %.1f", param[6]);
fprintf(printer,"\nCOEFFICIENT OF HISTERESIS(R) = %.1f", param[7]);
fprintf(printer,"\nTHE LOCAL ELASTIC MODULUS(IE) = %.2e",param[8]);
fprintf(printer,"\nTHE WORKING DIAMETER(diam) = %.2e",param[9]);
fprintf(printer,"\nTHE SEGMENT LENGTH(Z) = %.2e",param[10]);
fprintf(printer,"\nTHE DRAG COEFFICIENT(D) = %.4f", param[11]);
fprintf(printer,"\nTHE T-FORCE TO PROBABILITY SCALING(gain) = %.1f", param[12]);
fprintf(printer,"\nTHE ITERATION TIME INTERVAL(its) = %.5f",param[13]);
fputs("\n\n\tCURRENT MODE OF OPERATION", printer);
fprintf(printer,"\nTHE STATE OF THE MODEL(ACTIVE=1) = %.1f",param[14]);
fprintf(printer,"\nTHE BASE DRAG(1 highest; 0 lowest) = %.1f",param[15]);
fputs("\n\n\tTHIS IS PRINTED FROM THE C++ VERSION OF THE MODEL\n", printer);
fputc('\014', printer);

fclose(printer);
end  : ;
}//end of printout parameter function
//************************************
// find angles from slopes
double find_angle(double x, double y, double slope)
{
 double angle;
 if (slope==0.0) angle=0.0;
 else {
 if ( x>0 && y>0 )  angle=atan(slope);
 if ( x<0 && y>0 )  angle=atan(slope) + pi;
 if ( x<0 && y<0 )  angle=atan(slope) + pi;
 if ( x>0 && y<0 )  angle=atan(slope) + 2*pi;
      }
 return angle;
}
// function which graphs output on laser printer
void graph_output(void)
{
// local variables
int i,I,H,l;
double Z = param[10];
int LN = ceil(300000.0 * Z);
char choice;
char menu;
char date[10];
cout << "ENTER THE DATE; ";
cin >> date;
cout << "Assign an identifying code :"<< endl;
cout << "based on the simulated structure" << endl;
cout << "IDENTIFYING CODE:";
cin >> comm1;
printer =fopen ("prn", "w");
if (printer==0){
	 printf("printer was not located\n");
	 goto end;
	 }
fputs("\n\t ADVANCED MODEL in C version 7.1\n", printer);
fputs("\t   THE GEOMETRIC CLUTCH \n\n", printer);
fprintf(printer,"  OUTPUT DATED: %s    LENGTH IS: %d microns\n", date, LN);
fprintf(printer,"  IDENTIFYING CODE: %s\n", comm1);
fprintf(printer,"Plotted every %d iterations from %d to %d \n",Q,total-N,total);
printer =fopen ("prn", "wb");

fprintf(printer,"\33%0B");// graphic mode command
fprintf(printer,"IN");// initialize graphics Mode
fprintf(printer,"SP1;SC-700,700,-870,870;");
for(l=1; l<(PL+1); l++){
  fprintf(printer,"PA0,0");
for (i=1; i<31; i++){
	I = floor(G[l] [i]*(20/Z));
	H = floor(K[l] [i]*(20/Z));
	fprintf(printer,"PD%d,%d;",I,H);
	}
	if(l==1) fprintf(printer, "DT*,1;LB_s*;");
	if (l==PL) fprintf(printer,"LB_f*;");
	if (l==EX ) fprintf(printer,"LB_i*;");
	if ( EX>0 && l==(EX+1)) fprintf(printer,"LB_i+*;");
fprintf(printer,"PU"); // pen up
}
fprintf(printer,"SP0");
fprintf(printer,"\33%0A");// back to normal Mode
fprintf(printer,"\33E"); // reset & eject page
fclose(printer);
cout << "sent to printer" << endl;
getch();
end  : ;
}
//************************************
//************************************
// MAIN MODELING SUBROUTINE
void run_model(void)
{
// This begins the modelling portion of the Geometric Clutch program
// define local variables*********************
srand(time(NULL));
char C,c;
char Y,y;
int xval;
int yval;
int BRIDGE;
int COUNTER;
int i,j,k,n,l,A;
int RPTS=0;
const double DECAY_C=20000000000.0;
float R, BRVAL;
double TOTFORCE_R = .1;
double PANGL, LANGL, MAGN, XLEG, TORQUE, SUM_BP, SUM_BR, ADJUST,DC=0.0;
double AVEBCORR, TOTCOR, DRAGTORQUE, DYDX, ARM, VECTX, VECTY, PREVIOUS=0.0;
double HYPOT, DR, TF, PSI, SUMF, OLDDECAY,LONG_FORC_P, LONG_FORC_R, DECAY=0.0;
double TVFR, TVFP, OLFORCE, TOTFORCE_P, CHOICE, FEEDBK, NEWDECAY, PUSH=0.0;
double TRIGGER_P, TRIGGER_R, TRIGGEROOT_P, TRIGGEROOT_R, ADHESN_P, ADHESN_R=0.0;
double tfa[4] [31]={ {0},{0} };//t-force arrays 0 and 1 are long. and transverse comp.
//2 and 3 are P and R t-force at each segment
double adh[2] [32]={ {0} , {0} };// adhesion force storage
double forc[6] [31]={ {0} , {0} }; //longitudinal forces on the doublets
double dec[31]={0};// DECAY values
double osx[32]={0};// old position data from previous iteration
double osy[32]={0};
double drag[31]={0}; // dragtorque values
double torq[31]={0};// total torque passive and active
double equ[31]={0};// equilibrium curvatures
double adj[31]={0};// angle adjustments for feedback between iterations
PL =EX = 0;

// load parameters from global param array
reset :double coupl=param[1];
double TUG=param[2];
double SNAP=param[3];
double SWITCH_P=param[4];
double SWITCH_R=param[5];
double HISTER_P=param[6];
double HISTER_R=param[7];
double IE=param[8];
double DIAM=param[9];
double Z=param[10];
double D=param[11];
double GAIN=param[12];
double ITS=param[13];
double STATE=param[14];
double BASE=param[15];
//end of local variable list

repeat : A=1;
cls();
cout << "\t\t---THE GEOMETRIC CLUTCH MODEL---" ;
cout << "\n\n\n\n   Enter a number (100 to 1000) of iterations: ";
cin >> N;
cout << "Do you want an animation style screen display if yes type y:";
Y = getche();
y = toupper(Y);
if (y=='Y'){ Q = 1; goto anim;}

cout << "\n Enter a number(1 to 20) of iterations per plot: ";
cin >> Q;
anim : ;
cout << "\n To RUN press(ENTER)  To cancel type C: ";
C = getche();
c = toupper(C);
if (c=='C') goto last;
initgraph(&graphdriver, &graphmode, "c:\\tc\\bgi");
if (graphresult() < 0)
	exit(1);
setlinestyle (DASHED_LINE,0,NORM_WIDTH);
moveto (300,0);
lineto (300,460);
moveto (0,240);
lineto (640,240);
moveto (10,10);
outtext ("PRESS (ENTER) to START COMPUTING");
setlinestyle (SOLID_LINE,0,THICK_WIDTH);
if (y== 'Y') setwritemode (XOR_PUT);
getch();
for ( n=1; n<(N+1) ; n++)//MAIN ITERATION LOOP FOR MODELING
{

 if (n==1 || n%Q==0){ // Storage of position data for plotting
	PL+= 1;
	if (PL<49){
	for ( j=1; j<31; j++){
	G[PL] [j]=sx[j];
	K[PL] [j]=sy[j];
	}
     }
 }
if (n%Q==0 || n==1){
moveto (300,240);
if (n>1 && y=='Y'){
for ( i=1; i<31; i++)//ERASE OLD POSITION
	{
	xval= 300 + floor(osx[i] * 7.92/Z);
	yval= 240 - floor(osy[i] * 7.92/Z);
	lineto (xval,yval);
	}
     }
moveto (300,240);
for ( i=1; i<31; i++)//PLOT CURRENT POSITION
	{
	xval= 300 + floor(sx[i] * 7.92/Z);
	yval= 240 - floor(sy[i] * 7.92/Z);
	lineto (xval,yval);
	}
   }

// FIND T-FORCE
sang[1]=0;
for (j=1;j<31;j++){
	HYPOT =sqrt(((sang[j]*DIAM)*(sang[j]* DIAM))+pow(.0000033, 2));
	DR=HYPOT - .0000033;
	TF = SNAP * DR;
	if (sang[j]==0) sang[j]=.000001;
	PSI=fabs(atan(.0000033/(sang[j]*DIAM)));
	tfa[0] [j]=TF * cos(PSI) * (sang[j]/fabs(sang[j]));
	tfa[1] [j]=TF * sin(PSI);
	}
for ( l=1; l<31; l++) {//summation loop for passive force
	SUMF=0;
	for ( k=l; k<31; k++) {
		SUMF += tfa[0] [k];
		}
	forc[3] [l]=SUMF;
	}
if (STATE == 0) goto passive;
// FIND ACTIVE FORCES USING THE T-FORCE MECHANISM
for (l=1; l<31; l++) { // loop for bridge control
LONG_FORC_P = forc[3] [l] + forc[4] [l]*(1.0 - coupl) + coupl * forc[5] [l];
LONG_FORC_R = forc[3] [l] + forc[5] [l]*(1.0 - coupl) + coupl * forc[4] [l];
TVFP =tfa[1] [l] + Z*cur[l+1]*LONG_FORC_P;
TVFR =tfa[1] [l] + Z*cur[l+1]*LONG_FORC_R;
tfa[2] [l] = TVFP;
tfa[3] [l] = TVFR;

TRIGGEROOT_P = SWITCH_P + (TVFP * GAIN);
TRIGGEROOT_R = SWITCH_R + (TVFR * GAIN);

TRIGGER_R = TRIGGEROOT_R + adh[1] [l] +.5 * adh[1] [l-1]+ .5 * adh[1] [l+1];
if (TRIGGER_R<0.0) TRIGGER_R = 0.0;
TRIGGER_P = TRIGGEROOT_P + adh[0] [l] +.5 * adh[0] [l-1]+ .5 * adh[0] [l+1];
if (TRIGGER_P<0.0) TRIGGER_P = 0.0;

CHOICE = fabs(TOTFORCE_R);
if (fabs(TOTFORCE_P)>=fabs(TOTFORCE_R)) CHOICE = fabs(TOTFORCE_P);
if ( n%2==0 ) break;
OLFORCE = 1.0;
COUNTER = 0;
	while (fabs((CHOICE/OLFORCE)-1)>.1){
	COUNTER++;
	if (COUNTER>150) break;
	OLFORCE=CHOICE;
	TOTFORCE_P=0.0;
	TOTFORCE_R=0.0;
	  for (BRIDGE=1; BRIDGE<floor(4000000.0*Z); BRIDGE++){
	  R = rand();
	  BRVAL = R/RAND_MAX;
	  if (BRVAL<TRIGGER_P) TOTFORCE_P+= TUG;
	  R = rand();
	  BRVAL = R/RAND_MAX;
	  if (BRVAL<TRIGGER_R) TOTFORCE_R-= TUG;
	  }
	ADHESN_P = fabs(TOTFORCE_P) * ((1.0 - TRIGGEROOT_P) * HISTER_P);
	ADHESN_R = fabs(TOTFORCE_R) * ((1.0 - TRIGGEROOT_R) * HISTER_R);
	adh[0] [l] = ADHESN_P - coupl * ADHESN_R;
	adh[1] [l] = ADHESN_R - coupl * ADHESN_P;

TRIGGER_R = TRIGGEROOT_R + adh[1] [l] +.5 * adh[1] [l-1]+ .5 * adh[1] [l+1];
TRIGGER_P = TRIGGEROOT_P + adh[0] [l] +.5 * adh[0] [l-1]+ .5 * adh[0] [l+1];
	forc[1] [l] = TOTFORCE_P;
	forc[2] [l] = TOTFORCE_R;
	CHOICE = fabs(TOTFORCE_R);
	if (fabs(TOTFORCE_P)>=fabs(TOTFORCE_R)) CHOICE=fabs(TOTFORCE_P);
	if (OLFORCE == 0.0) OLFORCE=.0000001;
	if (CHOICE == 0) CHOICE = OLFORCE;
	}
}// end of bridge control loop
for ( i=1; i<31; i++){ // sum up forces from current position to tip
	SUM_BP=0;
	SUM_BR=0;
		for( j=i; j<31 ; j++){
		SUM_BP += forc[1] [j];
		SUM_BR += forc[2] [j];
		}
	forc[4] [i] = SUM_BP;
	forc[5] [i] = SUM_BR;
	} // END OF ACTIVE FORCE CALCULATION

passive  : ;
for ( l=1; l<31; l++) { //calculation of torques and equil curvatures
	forc[0] [l] = 2*forc[3] [l] + forc[4] [l] + forc[5] [l];
	torq[l]= forc[0] [l] * DIAM;
	equ[l]=-torq[l]/IE;
	OLDDECAY = dec[l];
	NEWDECAY = .6/(1+fabs(drag[l] * DECAY_C));
	if (n==1 && ROUTE==1) NEWDECAY=.06;
	if (n==2 && ROUTE==1) NEWDECAY=.04;
	if (n==3 && ROUTE==1) NEWDECAY=.02;
	DECAY = (OLDDECAY + NEWDECAY)/2;
	dec[l]= DECAY;
	adj[l]= (cur[l]-equ[l])*DECAY;
	}
for ( i=2; i<31; i++){ // Major feedback point
	cur[i]=cur[i] - adj[i];
	ang[i]=ang[i-1] + cur[i]*Z;
	sang[i]=ang[i] - ang[1];
	}
for ( j=1; j<31; j++){ // recalculate positions
	osx[j+1]=sx[j+1];
	osy[j+1]=sy[j+1];
	sx[j+1]=sx[j]+Z*cos(ang[j]);
	sy[j+1]=sy[j]+Z*sin(ang[j]);
	}                             // position reset loops
	// Drag routine for rotating base
dragloop : ;
if (BASE==1.0) goto maindrag;
AVEBCORR = TOTCOR = 0.0;
RPTS = 0;
drag[1] = .000001;
ADJUST = pow((Z/.0001), 3);
while ((fabs(drag[1]))>=(.000000014*ADJUST)){
 PREVIOUS = drag[1];
 RPTS+=1;
 if (RPTS>100) exit(1);
 l=1;
	DRAGTORQUE=0.0;
	DYDX=0.0;
	for (k = l+1; k<32; k++){
	ARM = sqrt(pow((sx[l]-sx[k]),2) + pow((sy[l]-sy[k]),2));
	VECTX = sx[k]-osx[k];
	VECTY = sy[k]-osy[k];
	if (VECTX==0) VECTX =.000000001;
	DYDX = VECTY/VECTX;
	PANGL = find_angle(VECTX, VECTY, DYDX);
	MAGN = sqrt(pow(VECTX, 2)+pow(VECTY,2));
	XLEG = sx[k]-sx[l];
	if (XLEG==0)  XLEG=.000000001;
	DYDX = (sy[k]-sy[l])/XLEG;
	LANGL= find_angle(XLEG, sy[k]-sy[l], DYDX);
	if (ARM==0) ARM = .000000001;
	TORQUE = ARM * D * Z *((MAGN * sin(PANGL-LANGL))/ITS);// individual contributions
	DRAGTORQUE += TORQUE; // running sum of individual contrib.
	   }
 drag[1]=DRAGTORQUE; // current value of drag at segment 1 the base

 FEEDBK = 6000/(D*ADJUST);
 PUSH = drag[1]*FEEDBK;
 AVEBCORR = -PUSH;
 TOTCOR+= AVEBCORR;
		for (j=1; j<31; j++) {
		 ang[j]+= AVEBCORR;
		 sx[j+1] = sx[j] + Z*cos(ang[j]);
		 sy[j+1] = sy[j] + Z*sin(ang[j]);
		 }
 if (RPTS>1 && (drag[1]*PREVIOUS)<0 ) drag[1]=0.0;
 }
	for (j=1; j<31; j++){
	ang[j]-=TOTCOR*BASE;
	sx[j+1]=sx[j]+Z*cos(ang[j]);
	sy[j+1]=sy[j]+Z*sin(ang[j]);
	}
maindrag : ;
for (l=1; l<31; l++){
	DRAGTORQUE=0.0;
	DYDX=0.0;
	for (k = l+1; k<32; k++){
	ARM = sqrt(pow((sx[l]-sx[k]),2) + pow((sy[l]-sy[k]),2));
	VECTX = sx[k]-osx[k];
	VECTY = sy[k]-osy[k];
	if (VECTX==0) VECTX =.000000001;
	DYDX = VECTY/VECTX;
	PANGL = find_angle(VECTX, VECTY, DYDX);
	MAGN = sqrt(pow(VECTX, 2)+pow(VECTY, 2));
	XLEG = sx[k]-sx[l];
	if (XLEG==0)  XLEG=.000000001;
	DYDX = (sy[k]-sy[l])/XLEG;
	LANGL= find_angle(XLEG, sy[k]-sy[l], DYDX);
	if (ARM==0) ARM = .000000001;
	DC = .5*D+(.5 * D * fabs(sin(PANGL-ang[k-1])));
	TORQUE = ARM * DC * Z *((MAGN * sin(PANGL-LANGL))/ITS);// individual contributions
	DRAGTORQUE += TORQUE; // running sum of individual contrib.
	   }
	drag[l]=DRAGTORQUE; // store them in the drag array
	}//end of drag loop

total+= 1;
}//END OF MAIN ITERATION LOOP
moveto(10, 400);
outtext("The run is finished:");
moveto(10, 420);
outtext("To continue press(ENTER)");
getch();
closegraph();
cout << "The model is now at iteration # " << total << endl;
cout << "Do you wish to continue?(if yes type(y), otherwise press(ENTER)): " << endl;
Y = getche();
y = toupper(Y);
if (y=='Y') {
	ROUTE = 2;
	cls();
	cout << "To continue computing press (ENTER)";
	PL=0;
	if (A==1) goto repeat;
	}
last : ;
ROUTE = 2;
}// END OF THE MODELING FUNCTION
// ********************************************************
// ************************* END OF PROGRAM ****************