#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define MAXYR  200
#define MAXEST 100
#define MSIZE  (MAXEST * (MAXEST + 1) / 2)

#define MAXSIM  600000
#define MAXSTP  500
#define CLCPAR  "CLC.PAR"
#define CLCDAT  "CLC.DAT"
#define CLCOUT  "CLC.OUT"

void RDPARS(void);
void PHOUT(double *CL, double RAWCL, int ILAST, int IY);
double CLIMIT(double *CATCH, double *SIGHT, double *FMATRX, int *ISYR, int *IZYR,
              double *ZMULT,double *POP, double *G);
double CONTRL(double D, double Y, double PLEVEL, double PSLOPE);
void DEVIAN(double *SS0, double *SS1, double *SS2, double *SS3, double *SIGHT, double *FMATRX, 
            int *ISYR, int *IZYR, double *ZMULT, double *POP, double *G);
double STKSIM(double RK, double R, double *POP, double *CATCH);
void SORT(double *ARRAY, double *ARRAY2, int N);
int  _getline(char *buff, size_t size, FILE *fp);

/* from Fortran COMMON VARIABLE */
struct MANPAR {
	double PPROB;  /* Probability level */
	double PYMAX;  /* Maximum value of the productivity parameter (Y) for integration */
	double PNYSTP; /* Number of step in integration over Y */
	double PKSTEP; /* Maximum relative step size in integration over K */
	double PDSTEP; /* Target step size for integration over depletion */
	double PNBSTP; /* Number of steps for integration over the bias parameter */
	double PBMIN;  /* Minimum multiplicative bias eg 0.5 means 50% downward bias */
	double PBMAX;  /* Maximum multiplicative bias eg 1.5 means 50% upward bias */
	double PSCALE; /* Raw deviance scaling factor S = 1/PSCALE**2 */
	double PHASET; /* Number of years without surveys before phaseout invoked */
	double PHASEP; /* Phaseout: anual reduction proportion */
	double PCYCLE; /* Maximum number of years before a new CLC = number of years a CLC is valid */
	double PLEVEL; /* Internal protection level */
	double PSLOPE; /* 'Slope' of catch control law */
} mp;

/* All years are scaled so 0 = a year prior to the 1st input data */
struct MANDAT {
	int ISTART;    /* Year of first input data (catch or abundance data) */
	int IYEAR;     /* Year for which to set first catch limit */
	int NS;        /* Number of non-zero estimates */
	int NZ;        /* Number of zero estimates */
	double RKLO;   /* Lower bound used in integration over K */
	double RKHI;   /* Upper bound used in integration over K */
} md;

int main(void)
{
/* ***** LOCAL VARIABLES */
/*
	CHARACTER STOCK*30, FORMT*50
	INTEGER IY, INITYR, I, IS, IYRCL, ILAST, N, IN, IOUT
*/
	double	CATCH[MAXYR + 1],  /* Catch array, indexed by year */
			SIGHT[MAXEST + 1], /* Abundance estimate array, indexed by estimate number */
			FMATRX[MSIZE + 1], /* Information matrix (H) of the log sightings estimates */
							 /* (excludeing zero estimates, stored as a lower triangle */
			ZMULT[MAXEST + 1], /* Poisson multiplier for Nth zero estimate */
			POP[MAXYR + 1],    /* Modelled population size in year I (set in STKSIM) */
			G[MAXYR + 1],      /* set & used in DEVIAN */
			RAWCL,           /* Nominal catch limit i.e. output from the Catch Limit Algorithm */
			CL;              /* The catch limit for the area considerd */
	int 	ISYR[MAXEST + 1],  /* Year of Nth abundance estimate SIGHT(N). N=0, NS-1 */
			IZYR[MAXEST + 1],  /* Year of Nth zero estimate. N=0, NZ-1 */
			POUT;            /* Parameter determing whether phaseout may be applied if */
							 /* necessary. Test for phaseout if POUT=1. (Phaseout is not */
							 /* applied to Medium of Large area nominal catch limits) */
	int 	INITYR,			 /* Year in which first premanagement catch taken */
			ILAST,			 /* Year of the most recent abundance estimate */
			IY, I, IS, IYRCL, N, J, N1;
	double	C, TOTALC;
	int		numlen, len, ret = 0;
	char /* STOCK[30], FORMT[50],*/ buff[300];

/* ***** Read in the data required for application of the CLA */

	/* DATA IN/2/, IOUT/3/ */
	FILE *IN = NULL, *IOUT = NULL;
	/* OPEN (IN, FILE = 'CLC.DAT') */
	if ((IN = fopen(CLCDAT, "rt")) == NULL) {
		fprintf(stderr, "Can't open " CLCDAT "\n");
		ret = 1;
		goto exit;
	}
	/* OPEN (IOUT, FILE = 'CLC.OUT') */
	if ((IOUT = fopen(CLCOUT, "w")) == NULL) {
		fprintf(stderr, "Can't open " CLCOUT "\n");
		ret = 1;
		goto exit;
	}

	RDPARS();

	/* READ (IN, '(A30/)') STOCK */
	_getline(buff, sizeof(buff), IN);
	/* WRITE (IOUT, '(A30/)') STOCK */
	fprintf(IOUT, "%s by C\n\n", buff);
	fprintf(stdout, "%s by C\n", buff);
	_getline(buff, sizeof(buff), IN); /* Discard empty lines */

/* 	Read the year of the first catch, the first year for which catch */
/* limits are to be set & the phaseout option */

	/* READ (IN, '(T30,I4)') INITYR */
	_getline(buff, sizeof(buff), IN);
	numlen = 4; len = strlen(buff);
	INITYR = atoi(&buff[len-numlen]);

	/* WRITE (IOUT, '(A,T30,I4)') 'Year of first input data', INITYR */
	fprintf(IOUT,   "%.*s %4d\n", len - numlen, buff, INITYR);
	fprintf(stdout, "%.*s %4d\n", len - numlen, buff, INITYR);

	/* READ (IN, '(T30,I4)') IYRCL */
	_getline(buff, sizeof(buff), IN);
	numlen = 4; len = strlen(buff);
	IYRCL = atoi(&buff[len-numlen]);

	/* WRITE (IOUT, '(A, T30, I4)') 'Year of first catch limit', IYRCL */
	fprintf(IOUT,   "%.*s %4d\n", len - numlen, buff, IYRCL);
	fprintf(stdout, "%.*s %4d\n", len - numlen, buff, IYRCL);

	/* READ (IN, '(T30, I4)') POUT */
	_getline(buff, sizeof(buff), IN);
	numlen = 4; len = strlen(buff);
	POUT = atoi(&buff[len-numlen]);
	if (POUT == 1) {
		/* WRITE (IOUT, '(A, T30,A4)') 'Apply phaseout if necessary', 'YES' */
		fprintf(IOUT,   "Apply phaseout if necessary    YES\n");
		fprintf(stdout, "Apply phaseout if necessary    YES\n");
	} else {
		/* WRITE (IOUT, '(A, T30,A4)') 'Apply phaseout', 'No' */
		fprintf(IOUT,   "Apply phaseout                  No\n");
		fprintf(stdout, "Apply phaseout                  No\n");
	}

/* Re-scale IYRCL such that 0 is the year prior to the first input data */

	md.ISTART = 0;
	md.IYEAR = IYRCL - INITYR;
	/* IF (IYEAR .LE. 0 .OR. IYEAR .GT. MAXYR) STOP 'INVALID YEAR' */
	if ((md.IYEAR <= 0) || (md.IYEAR > MAXYR)) {
		fprintf(stderr, "INVALID YEAR");
		ret = 1;
		goto exit;
	}

/* Initialize the catch array */

	/* DO 10 I = 0, MAXYR */
	for (I = 0; I <= MAXYR; I++) {
		CATCH[I] = 0.0;
	} /* 10   CONTINUE */

/* Read in the catch data, scaling each year to the initial year */
	/* READ (IN, '(// A)') FORMT (I4,F8.0) */
	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
	_getline(buff, sizeof(buff), IN); /* Discard empty lines */

	/* WRITE (IOUT, '(/A)') 'Historic catches:' */
	fprintf(IOUT, "\nHistoric catches:\n");
	fprintf(stdout, "\nHistoric catches:\n");
	TOTALC = 0;
	/* DO 20 I = 0, MAXYR */
	for (I = 0; I <= MAXYR; I++) {
		/* READ (IN, FORMT) IY, C */ 
		_getline(buff, sizeof(buff), IN);
		sscanf(buff,"%4d%7lf", &IY, &C);
		/* IF (IY .LT. 0) GO TO 25 */
		if (IY < 0) break; /* EOD */

		/* WRITE (IOUT, FORMT) IY, C */
		fprintf(stdout, "%4d%6.f\n", IY, C);
		fprintf(IOUT, "%4d%6.f\n", IY, C);
		IY = IY - INITYR;
		/* IF (IY .LT. 0 .OR. IY .GE. IYEAR) STOP */
		if (IY < 0 || IY > md.IYEAR) {
			fprintf(stderr, " ERROR: IY is out of range(%d)\n", IY);
			ret = 1;
			goto exit;
		}
		CATCH[IY] = C;
		TOTALC = TOTALC + C;
	} /* 20   CONTINUE */

/* 25   IF (TOTALC .LT. 0) STOP ' ERROR: No historic catch input' */
	if (TOTALC < 0) {
		fprintf(stderr, " ERROR: No historic catch input\n");
		ret = 1;
		goto exit;
	}

/* Read in the non-zero sightings estimates and information matrix */

	/* READ (IN, '(//T30, I4 / A)') NS, FORMT (I4,F10.0,10F8.0) */
	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
	_getline(buff, sizeof(buff), IN);
	numlen = 4; len = strlen(buff);
	md.NS = atoi(&buff[len-numlen]);

	/* WRITE (IOUT, '(/A)') 'Abundance estimates:' */
	fprintf(IOUT,   "\nAbundance estimates:\n");
	fprintf(stdout, "\nAbundance estimates:\n");
	/* IF (NS .GT. MAXEST) STOP 'ERROR: Abundance year out of range' */
	if (md.NS > MAXEST) {
		fprintf(stderr, "ERROR: Abundance year out of range\n");
		ret = 1;
		goto exit;
	}
	_getline(buff, sizeof(buff), IN);
	/* DO 30 N=0, NS-1 */
	for (N = 0; N <= md.NS - 1; N++) {
		N1 = (N * (N + 1)) / 2;
		/* READ (IN, FORMT) ISYR(N), SIGHT(N), (FMATRX(N1 + J), J = 0, N) */
		_getline(buff, sizeof(buff), IN);
		sscanf(buff, "%4d%10lf", &ISYR[N], &SIGHT[N]);
		fprintf(IOUT,   "%4d %8.f", ISYR[N], SIGHT[N]);
		fprintf(stdout, "%4d %8.f", ISYR[N], SIGHT[N]);
		len = strlen(buff);
		for (J = 0; J < N + 1; J++) {
			if (len > J * 7 + 14) {
				sscanf(&buff[0] + J * 7 + 14,"%7lf", &FMATRX[N1 + J]);
			} else {
				FMATRX[N1 + J] = 0.0;
			}
			fprintf(IOUT,   "%7.f", FMATRX[N1 + J]);
			fprintf(stdout, "%7.f", FMATRX[N1 + J]);
		}
		fprintf(IOUT,   "\n");
		fprintf(stdout, "\n");
		/* WRITE (IOUT, FORMT) ISYR(N), SIGHT(N), (FMATRX (N1 + J), J = 0, N) */
		ISYR[N] = ISYR[N] - INITYR;
		if (ISYR[N] < 0 || ISYR[N] >= md.IYEAR) {
			/*STOP 'ERROR: Sight year out of range'*/
			fprintf(stderr, "ERROR: Sight year out of range\n");
			ret = 1;
			goto exit;
		}
		if (SIGHT[N] <= 0.0) {
			/* STOP ' ERROR: Estimate not positive' */
			fprintf(stderr, "ERROR: Estimate not positive\n");
			ret = 1;
			goto exit;
		}
	} /* 30   CONTINUE */

/* Read in the Poisson multipliers for any zero sightings estimates */

	/* READ (IN, '(/T30, I3, /A)') NZ, FORMAT */
	_getline(buff, sizeof(buff), IN);
	_getline(buff, sizeof(buff), IN);
	numlen = 3; len = strlen(buff);
	md.NZ = atoi(&buff[len-numlen]);
	if (md.NZ > 0) {
		/* WRITE (IOUT, '(/A)') 'Zero abundance estimates:' */
		fprintf(IOUT,   "\nZero abundance estimates:\n");
		fprintf(stdout, "\nZero abundance estimates:\n");
	}
	if (md.NZ > MAXEST) {
		/* STOP ' ERROR: Zero estimate array too small' */
		fprintf(stderr, "ERROR: Zero estimate array too small\n");
		ret = 1;
		goto exit;
	}
	_getline(buff, sizeof(buff), IN);
	/* DO 40 N = 0, NZ - 1 */
	for (N = 0; N <= md.NZ - 1; N++) {
		N1 = (N * (N + 1)) / 2;
		/*READ (IN, FORMT) IZYR(N), ZMULT(N)*/
		_getline(buff, sizeof(buff), IN);
		sscanf(buff, "%4d%10lf", &IZYR[N], &ZMULT[N]);
		/* WRITE (IOUT, FORMT) IZYR(N), ZMULT(N) */
		fprintf(IOUT,   "%4d%7.f\n", IZYR[N], ZMULT[N]);
		fprintf(stdout, "%4d%7.f\n", IZYR[N], ZMULT[N]);
		IZYR[N] = IZYR[N] - INITYR;
		if (IZYR[N] < 0 || IZYR[N] >= md.IYEAR) {
			/* STOP ' Sight year out of range' */
			fprintf(stderr, "Sight year out of range\n");
			ret = 1;
			goto exit;
		}
                                    
		if (ZMULT[N] <= 0.0) {
			/* STOP ' ERROR: Multiplier not positive' */
			fprintf(stderr, "ERROR: Multiplier not positive\n");
			ret = 1;
			goto exit;
		}
	} /* 40   CONTINUE */
	/*WRITE (IOUT, '()')*/
	fprintf(IOUT,"\n");
	fprintf(stdout,"\n");
/* Bound the range for the integration over K */

	md.RKHI = 1.0E7;
	md.RKLO = 0.0;
	IS = 0; /* orignal not initialized bug? */
/* ****** */
/* ****** Run the CLA to obtain the nominal catch limit */
/* ****** */
	RAWCL = CLIMIT(CATCH, SIGHT, FMATRX, ISYR, IZYR, ZMULT, POP, G);
/* Set the catch limits for PCYCLE years. If the catch limit may be */
/* subject to phaseout, call POUT to apply the phaseout rule. */

/* First set ILAST = year of the most recent abundance estimate */

	if (md.NS > 0) ILAST = ISYR[md.NS - 1];
	/* IF (NZ .GT. 0) ILAST = MAX(IS, IZYR(md.NZ - 1)) */
	if (md.NZ > 0) ILAST = (IS > IZYR[md.NZ - 1])? IS: IZYR[md.NZ - 1];

	/* DO 100 IY = IYEAR, IYEAR + mp.PCYCLE - f1 */
	for (IY = md.IYEAR; IY <= md.IYEAR + mp.PCYCLE - 1; IY++) {
		if (POUT == 1) {
			PHOUT(&CL, RAWCL, ILAST, IY);
		} else {
			CL = RAWCL;
		}
		/* WRITE (IOUT, '(A6, I5, A17, I6)') 'Year:', IY + INITYR,'Catch limit:', NINT(CL) */
		fprintf(IOUT,   "Year: %4d     Catch limit:%6d\n", IY + INITYR, (int)(CL + 0.5));
		fprintf(stdout, "Year: %4d     Catch limit:%6d\n", IY + INITYR, (int)(CL + 0.5));

	} /* 100  CONTINUE */

exit:
	if (IN != NULL)   fclose(IN);
	if (IOUT != NULL) fclose(IOUT);
	return ret;
	/*	STOP */
	/*	END  */
}

/*
C **************************************************************************************
C
C	CLC version 6
C
C **************************************************************************************
C
C	31 January 1994
C
C **************************************************************************************
*/
/*
	SUBROUTINE RDPARS
	Read the file of input parameters
*/
void RDPARS(void)
{
	double *p[] = {
		&mp.PPROB, &mp.PYMAX, &mp.PNYSTP,&mp.PKSTEP, &mp.PDSTEP, &mp.PBMIN, &mp.PBMAX,
		&mp.PNBSTP, &mp.PSCALE, &mp.PHASET, &mp.PHASEP, &mp.PCYCLE, &mp.PLEVEL, &mp.PSLOPE
	};
	FILE *IFILE;
	char buff[80];
	int count;
	
	/* Open the input file (only read once) */
	IFILE = fopen(CLCPAR, "rt");
	if (IFILE == NULL) {
		fprintf(stderr,"Can't open " CLCPAR "\n");
		exit(1);
	}
	count = 0;
	while (_getline(buff, sizeof(buff), IFILE) != 0) {
		/* READ (IFILE, '(T30, F10.0)') PPROB, PYMAX, PNYSTP, PKSTEP, PDSTEP, PBMIN, PBMAX,
								 PNBSTP, PSCALE, PHASET, PHASEP, PCYCLE, PLEVEL, PSLOPE */
		sscanf(buff,"%*28c%lf", p[count]);
		count++;
	}
	/* CLOSE (IFILE) */
	fclose(IFILE);
	return;
	/* RETURN
	   END */
}

/*
	SUBROUTINE PHOUT (CL, RAWCL, ILAST, IY)
*/
void PHOUT(double *CL, double RAWCL, int ILAST, int IY)
{
/*
	PHASET Number of years without surveys before phaseout invoked
	PHASEP Phaseout annual reduction proportion
	Phaseout: Reduce catch limit if there is no survey data in the
	last PHASET years
*/
	if (IY >= ILAST + mp.PHASET) {
		*CL = RAWCL * (1.0 - mp.PHASEP * (IY - ILAST - mp.PHASET));
		if (*CL < 0.0) *CL = 0.0;
	} else {
		*CL = RAWCL;
	}
	/* RETURN
	   END */
}

/*
	REAL FUNCTION CLIMIT (CATCH,SIGHT,FMATRX,ISYR,IZYR,ZMULT,POP,G)
	Run the CLA to obtain to nominal catch limit
*/
double CLIMIT(double *CATCH, double *SIGHT, double *FMATRX, int *ISYR, int *IZYR, double *ZMULT,double *POP, double *G)
{
/*
	Local variables: 

	REAL PRES(0:MAXSIM), QRES(0:MAXSIM), SS0, SS1, SS2, SS3
	REAL SF, Y(0:MAXSTP), B(0:MAXSTP), RLGB(0:MAXSTP)
	INTEGER NB, NR
	NB	Number of bias steps
	NR	Number of steps for Y (the productivity parameter)
	SF	Deviance scale factor (SF = .5 / PSCALE**2)
*/
	double SF, SS0, SS1, SS2, SS3;
	double Y[MAXSTP + 1], B[MAXSTP + 1], RLGB[MAXSTP + 1];
	double *PRES = NULL, *QRES = NULL;

	int NB, NR;
	double BINC, YINC, DK, R, RK, D, DD, P, Q, PTOT;
	int I, J, N, N2;

	double _CLIMIT = 0.0;
	if (md.NS <= 0) return _CLIMIT;

	PRES = (double *) malloc(sizeof(double) * (MAXSIM + 1));
	if (PRES == NULL) {
		fprintf(stderr, "Can't Allocate memory (PRES)\n");
		goto err_exit;
	}
	memset(PRES, 0, sizeof(double) * (MAXSTP + 1));

	QRES = (double *) malloc(sizeof(double) * (MAXSIM + 1));
	if (QRES == NULL) {
		fprintf(stderr, "Can't Allocate memory (QRES)\n");
		goto err_exit;
	}
	memset(QRES, 0, sizeof(double) * (MAXSTP + 1));

	/* Set deviance scale factor S = 1/PSCALE**2 */
	SF = 0.5 / (mp.PSCALE * mp.PSCALE);

	/* Check the sizes of the Y and B arrays are large enough */
	if (mp.PNBSTP > MAXSTP || mp.PNYSTP > MAXSTP) {
		fprintf(stderr, "Y &/or b array sizes too small\n");
		goto err_exit;
	}
	/* Set sightings bias step sizes & their log values. BINC = Bias increment */
	NB = (int)mp.PNBSTP;
	BINC = (mp.PBMAX - mp.PBMIN) / NB;
	/* DO 10 I = 0, NB - 1 */
	for (I = 0; I <= NB - 1; I++) {
        *(B + I) = mp.PBMIN + (I + 0.5) * BINC;
        *(RLGB + I) = -log(*(B + I));
	} /* 10   CONTINUE */

	/* Set productivity parameter step sizes (midpoints) */
	NR = (int)mp.PNYSTP;
	YINC = mp.PYMAX / NR;
	/* DO 20 I = 0, NR - 1 */
	for (I = 0; I <= NR - 1; I++) {
		*(Y + I) = (I + 0.5) * YINC;
	}/* 20   CONTINUE */

	PTOT = 0;
	N = 0;
	/* DO 50 I = 0, NR - 1 */
	for (I = 0; I <= NR - 1; I++) {
		/* Set R from the productivity parameter, Y */
		R = 1.4184 * *(Y + I);

		/* Step size for K */
        DK = mp.PKSTEP;
        D = 1.0;
        RK = md.RKHI;

		/* Use function STKSIM to set up the Nth population trajectory */
		/* i.e. set up the pop array */

		while (1) {
		/* 30     IF (RK .LE. RKLO .OR. STKSIM (RK, R, POP, CATCH) .LE. 0.) GOTO 40 */
			if ((RK <= md.RKLO) || (STKSIM(RK, R, POP, CATCH) <= 0.0)) break;
			/* IF (N .GE. MAXSIM) STOP 'ERROR: TOO MANY SIMULATIONS' */
			if (N >= MAXSIM) {
				fprintf(stderr, "ERROR: TOO MANY SIMULATIONS");
				goto err_exit;
			}
			/* How much depletion covered? */
			DD = D - POP[md.IYEAR] / RK;
			D = POP[md.IYEAR] / RK;
			P = 0.0;

			if (DD > 0.0) {
				/* Compute the internal catch limit corresponding to D and Y(I) */
				QRES[N] = CONTRL(D, Y[I], mp.PLEVEL, mp.PSLOPE) * POP[md.IYEAR];

				/* Calculate deviance */
            	DEVIAN(&SS0, &SS1, &SS2, &SS3, SIGHT, FMATRX, ISYR,IZYR, ZMULT, POP, G);

				/* Scale likelihood and integrate over values for the bias parameter */
				/* DO 35 J = 0, NB - 1 */
				for (J = 0; J <= NB - 1; J++) {
					P = P + exp(-SF * (SS0 + *(RLGB + J) * (SS1 + *(RLGB + J) * SS2) + SS3 * *(B + J)));
				} /* 35         CONTINUE */

				/* Calculate the weight for this point (& total area under likelihood) */
				PRES[N] = P * DD;
				PTOT = PTOT + PRES[N];

				/* Update counter */
				N = N + 1;

				/* Find the next K */
				DK = DK * mp.PDSTEP / DD;
				/* IF (DK .GT. PKSTEP) DK = PKSTEP */
				if (DK > mp.PKSTEP) DK = mp.PKSTEP;

			} else {
				/* IF DD = 0 change the step size only */
				DK = mp.PKSTEP;
			}
			/* Set the new value of K */
			RK = RK / (1.0 + DK);
		} /* GOTO 30 */ /* 40     CONTINUE */
	} /* 50   CONTINUE */
	/* IF (PTOT .LE. 0.) STOP 'ERROR: PROB INTEGRATES TO ZERO' */
	if (PTOT < 0.0) {
		fprintf(stderr,"ERROR: PROB INTEGRATES TO ZERO");
		goto err_exit;
	}

	/* Sort the QRES and PRES arrays in ascending order of QRES */
	N2 = N;
	SORT(QRES, PRES, N2);

	/* Normalize the relative likelihoods */
	/* DO 60 I = 0, N - 1 */
	for (I = 0; I <= N - 1; I++) {
		PRES[I] = PRES[I] / PTOT;
	} /* 60   CONTINUE */

	/* Extract the desired probability level: the nominal catch limit (NCL) */
	/* is the lower PROB% of the distribution. */
	/* First calculate PRES(I), the probability that the NCL is between */
	/* QRES(I) & QRES(I + 1). */

	P = 0;
	/* DO 70 I = 0, N - 1 */
	for (I = 0; I <= N - 1; I++) {
		P = P + PRES[I];
		/* IF (P .GT. PPROB) GOTO 80 */
		if (P > mp.PPROB) break;
	} /* 70   CONTINUE */

	/* Interpolate to set the nominal catch limit */
	/* 80   IF (I .GE. N - 1) THEN */
	if  (I >= N - 1) {
        Q = QRES[N - 1];
	} else {
        Q = (QRES[I + 1] * (mp.PPROB - P + PRES[I]) + QRES[I] * (P - mp.PPROB)) / PRES[I];
	}

	_CLIMIT = Q;

	if (QRES != NULL) free(QRES);
	if (PRES != NULL) free(PRES);
/*	if (Y != NULL) free(Y);
	if (B != NULL) free(B);
	if (RLGB != NULL) free(RLGB);*/
	return _CLIMIT;

err_exit:
	if (QRES != NULL) free(QRES);
	if (PRES != NULL) free(PRES);
/*	if (Y != NULL) free(Y);
	if (B != NULL) free(B);
	if (RLGB != NULL) free(RLGB);*/
	exit(1);
	return 0.0; /* not reached */
	/* RETURN */
	/* END */
}

/*
	FUNCTION CONTRL (D, Y, PLEVEL, PSLOPE)
	Catch control law
*/
double CONTRL(double D, double Y, double PLEVEL, double PSLOPE)
{
	double ret;
	if (D < PLEVEL) {
		ret = 0.0;
	} else if (D < 1.0) {
		ret = PSLOPE * Y * (D - PLEVEL);
	} else {
		ret = PSLOPE * Y * (1.0 - PLEVEL);
	}
	return ret;
	/* END */
}

/*
	SUBROUTINE DEVIAN (SS0, SS1, SS2, SS3, SIGHT, FMATRX, ISYR,
	Calculate deviance (-2 log likelihood) in terms of coefficients for 
	the bias and log bias parameters
*/
void DEVIAN(double *SS0, double *SS1, double *SS2, double *SS3, double *SIGHT,
            double *FMATRX, int *ISYR, int *IZYR, double *ZMULT, double *POP, double *G)
{
	/* REAL SS0,SS1,SS2,SIGHT(0:*),FMATRX(0:*),ZMULT(0:*),POP(0:*),G(0:*) */
	/* INTEGER ISYR(0:*), IZYR(0:*) */
	int N, J, K;

	*SS0 = 0.0;
	*SS1 = 0.0;
	*SS2 = 0.0;
	*SS3 = 0.0;

	/* DO 100 N = 0, NS - 1 */
	for (N = 0; N <= md.NS - 1; N++) { 

		*(G + N) = log(*(SIGHT + N) / *(POP+*(ISYR + N)));
		K = N * (N + 1) / 2;

		/* DO 10 J = 0, N - 1 */
		for (J = 0; J <= N - 1; J++) {
			/*1st add non diagonal contributions (which are doubled up)*/
			*SS0 = *SS0 + 2.0 * *(G + J) * *(G + N) * *(FMATRX + K);
			*SS1 = *SS1 + 2.0 * (*(G + J) + *(G + N)) * *(FMATRX + K);
			*SS2 = *SS2 + *(FMATRX + K) + *(FMATRX + K);
			K = K + 1;
		} /* 10     CONTINUE*/

		/* Now add diagnoal contribution */
        *SS0 = *SS0 + *(G + N) * *(G + N) * *(FMATRX + K);
        *SS1 = *SS1 + 2.0 * *(G + N) * *(FMATRX + K);
        *SS2 = *SS2 + *(FMATRX + K);

	} /* 100  CONTINUE */

	/* Now do the zero estimates */
	/* DO 200 N = 0, NZ - 1 */
	for (N =0; N <= md.NZ - 1; N++) {
		*SS3 = *SS3 + 2.0 * *(POP + *(IZYR + N)) / *(ZMULT + N);
 	} /* 200  CONTINUE */

	/* RETURN */
	/* END    */
	return;
}
/*
	FUNCTION STKSIM (RK, R, POP, CATCH)
	Calculate the stock trajectory with parameter RK and R: return
	the current stock size
	RK = notional carrying capacity; R = productivity parameter * 1.4184 
*/
double STKSIM(double RK, double R, double *POP, double *CATCH)
{
	/*INTEGER ISTART, IYEAR, NS, NZ, RKLO, RKHI */
	/* REAL CATCH(0:*), POP(0:*) */
	int I;
	double D;
	*(POP+md.ISTART) = RK;
	/*DO 10 I = md.ISTART + 1, md.IYEAR*/
	for (I = md.ISTART + 1; I <= md.IYEAR; I++) {
		D = *(POP + I - 1) / RK;
		*(POP + I) = *(POP + I - 1) * (1.0 + R * (1.0 - D * D)) - *(CATCH + I -1);
		if (*(POP + I) <= 0) {	/*GOTO 20*/
			return 0.0;
			/* 20 STKSIM = 0. */
		}
	} /* 10   CONTINUE */
	/* STKSIM = POP(md.IYEAR) */
	return *(POP + md.IYEAR);
	/* RETURN */
    /* END */
}
/* 
	SUBROUTINE SORT (ARRAY, ARRAY2, N)
	SORT sorts a pair of arrays in ascending order of the first array
   (using the heapsort algorithm)
*/
void SORT(double *ARRAY, double *ARRAY2, int N)
{
	double /*ARRAY(0:*), ARRAY2(0:*),*/ TEMP, TEMP2;
	int  /*N,*/ K, IR, I, J;
	if (N < 2) return;
	K = N / 2;
	IR = N - 1;
	while (1) { /*10*/
		if (K != 0) {
			K = K - 1;
			TEMP = *(ARRAY + K);
			TEMP2 = *(ARRAY2 + K);
		} else {
			TEMP = *(ARRAY + IR);
			TEMP2 = *(ARRAY2 + IR);
			*(ARRAY + IR) = *ARRAY;
			*(ARRAY2 + IR) = *ARRAY2;
			IR = IR - 1;
			if (IR == 0) {
				*ARRAY = TEMP;
				*ARRAY2 = TEMP2;
				return;
			}
		}
		I = K;
		J = K + K + 1;
		while (1) {/*20*/
			if (J > IR) break;
			if ((J < IR) && (*(ARRAY + J) < *(ARRAY+ J + 1))) J = J + 1;
			if (TEMP < *(ARRAY + J)) {
				*(ARRAY + I) = *(ARRAY + J);
				*(ARRAY2 + I) = *(ARRAY2 + J);
				I = J;
				J = J + I + 1;
			} else {
				J = IR + 1;
			}
		} /*GOTO 20*/
		*(ARRAY + I) = TEMP;
		*(ARRAY2 + I) = TEMP2;
	}/*GOTO 10*/
	/*END*/
}
int _getline(char *buff, size_t size, FILE *fp){
	char *ret;
	ret = fgets(buff, size, fp);
	if (ret != NULL) {
		if (buff[strlen(buff) - 1] == '\n')
			buff[strlen(buff) - 1] = 0;
		if (buff[strlen(buff) - 1] == '\r')
			buff[strlen(buff) - 1] = 0;
	}
	return (ret != NULL);
}