fork download
copy 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include <math.h>
  5.  
  6. #define MAXYR  200
  7. #define MAXEST 100
  8. #define MSIZE  (MAXEST * (MAXEST + 1) / 2)
  9.  
  10. #define MAXSIM  600000
  11. #define MAXSTP  500
  12.  
  13. #define IN stdin
  14. #define IOUT stdout
  15.  
  16. void PHOUT(double *CL, double RAWCL, int ILAST, int IY);
  17. double CLIMIT(double *CATCH, double *SIGHT, double *FMATRX, int *ISYR, int *IZYR,
  18.               double *ZMULT,double *POP, double *G);
  19. double CONTRL(double D, double Y, double PLEVEL, double PSLOPE);
  20. void DEVIAN(double *SS0, double *SS1, double *SS2, double *SS3, double *SIGHT, double *FMATRX, 
  21.             int *ISYR, int *IZYR, double *ZMULT, double *POP, double *G);
  22. double STKSIM(double RK, double R, double *POP, double *CATCH);
  23. void SORT(double *ARRAY, double *ARRAY2, int N);
  24. int  _getline(char *buff, size_t size, FILE *fp);
  25.  
  26. /* from Fortran COMMON VARIABLE */
  27. struct MANPAR {
  28. 	double PPROB;  /* Probability level */
  29. 	double PYMAX;  /* Maximum value of the productivity parameter (Y) for integration */
  30. 	double PNYSTP; /* Number of step in integration over Y */
  31. 	double PKSTEP; /* Maximum relative step size in integration over K */
  32. 	double PDSTEP; /* Target step size for integration over depletion */
  33. 	double PNBSTP; /* Number of steps for integration over the bias parameter */
  34. 	double PBMIN;  /* Minimum multiplicative bias eg 0.5 means 50% downward bias */
  35. 	double PBMAX;  /* Maximum multiplicative bias eg 1.5 means 50% upward bias */
  36. 	double PSCALE; /* Raw deviance scaling factor S = 1/PSCALE**2 */
  37. 	double PHASET; /* Number of years without surveys before phaseout invoked */
  38. 	double PHASEP; /* Phaseout: anual reduction proportion */
  39. 	double PCYCLE; /* Maximum number of years before a new CLC=number of years a CLC is valid */
  40. 	double PLEVEL; /* Internal protection level */
  41. 	double PSLOPE; /* 'Slope' of catch control law */
  42. } mp = {
  43. 	0.4102, 0.05,  200.0,  0.05,  0.01,
  44. 	100.0,  0.0,   1.6667, 4.0,   8.0,
  45. 	0.2,    5.0,   0.54,   3.0
  46. };
  47. /* All years are scaled so 0 = a year prior to the 1st input data */
  48. struct MANDAT {
  49. 	int ISTART;    /* Year of first input data (catch or abundance data) */
  50. 	int IYEAR;     /* Year for which to set first catch limit */
  51. 	int NS;        /* Number of non-zero estimates */
  52. 	int NZ;        /* Number of zero estimates */
  53. 	double RKLO;   /* Lower bound used in integration over K */
  54. 	double RKHI;   /* Upper bound used in integration over K */
  55. } md;
  56.  
  57. int main(void)
  58. {
  59. /* ***** LOCAL VARIABLES */
  60. /*
  61. 	CHARACTER STOCK*30, FORMT*50
  62. 	INTEGER IY, INITYR, I, IS, IYRCL, ILAST, N, IN, IOUT
  63. */
  64. 	double	CATCH[MAXYR + 1],  /* Catch array, indexed by year */
  65. 			SIGHT[MAXEST + 1], /* Abundance estimate array, indexed by estimate number */
  66. 			FMATRX[MSIZE + 1], /* Information matrix (H) of the log sightings estimates */
  67. 							 /* (excludeing zero estimates, stored as a lower triangle */
  68. 			ZMULT[MAXEST + 1], /* Poisson multiplier for Nth zero estimate */
  69. 			POP[MAXYR + 1],    /* Modelled population size in year I (set in STKSIM) */
  70. 			G[MAXYR + 1],      /* set & used in DEVIAN */
  71. 			RAWCL,           /* Nominal catch limit i.e. output from the Catch Limit Algorithm */
  72. 			CL;              /* The catch limit for the area considerd */
  73. 	int 	ISYR[MAXEST + 1],  /* Year of Nth abundance estimate SIGHT(N). N=0, NS-1 */
  74. 			IZYR[MAXEST + 1],  /* Year of Nth zero estimate. N=0, NZ-1 */
  75. 			POUT;            /* Parameter determing whether phaseout may be applied if */
  76. 							 /* necessary. Test for phaseout if POUT=1. (Phaseout is not */
  77. 							 /* applied to Medium of Large area nominal catch limits) */
  78. 	int 	INITYR,			 /* Year in which first premanagement catch taken */
  79. 			ILAST,			 /* Year of the most recent abundance estimate */
  80. 			IY, I, IS, IYRCL, N, J, N1;
  81. 	double	C, TOTALC;
  82. 	int		numlen, len, ret = 0;
  83. 	char /* STOCK[30], FORMT[50],*/ buff[300];
  84.  
  85. /* ***** Read in the data required for application of the CLA */
  86.  
  87. 	_getline(buff, sizeof(buff), IN);
  88. 	fprintf(IOUT, "%s by C\n\n", buff);
  89. 	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
  90.  
  91. /* 	Read the year of the first catch, the first year for which catch */
  92. /* limits are to be set & the phaseout option */
  93.  
  94. 	_getline(buff, sizeof(buff), IN);
  95. 	numlen = 4; len = strlen(buff);
  96. 	INITYR = atoi(&buff[len-numlen]);
  97.  
  98. 	fprintf(IOUT,   "%.*s %4d\n", len - numlen, buff, INITYR);
  99.  
  100. 	_getline(buff, sizeof(buff), IN);
  101. 	numlen = 4; len = strlen(buff);
  102. 	IYRCL = atoi(&buff[len-numlen]);
  103.  
  104. 	fprintf(IOUT,   "%.*s %4d\n", len - numlen, buff, IYRCL);
  105.  
  106. 	_getline(buff, sizeof(buff), IN);
  107. 	numlen = 4; len = strlen(buff);
  108. 	POUT = atoi(&buff[len-numlen]);
  109. 	if (POUT == 1) {
  110. 		fprintf(IOUT,   "Apply phaseout if necessary    YES\n");
  111. 	} else {
  112. 		fprintf(IOUT,   "Apply phaseout                  No\n");
  113. 	}
  114.  
  115. /* Re-scale IYRCL such that 0 is the year prior to the first input data */
  116.  
  117. 	md.ISTART = 0;
  118. 	md.IYEAR = IYRCL - INITYR;
  119. 	if ((md.IYEAR <= 0) || (md.IYEAR > MAXYR)) {
  120. 		fprintf(stderr, "INVALID YEAR");
  121. 		ret = 1;
  122. 		goto exit;
  123. 	}
  124.  
  125. /* Initialize the catch array */
  126.  
  127. 	for (I = 0; I <= MAXYR; I++) {
  128. 		CATCH[I] = 0.0;
  129. 	}
  130.  
  131. /* Read in the catch data, scaling each year to the initial year */
  132. 	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
  133. 	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
  134. 	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
  135.  
  136. 	fprintf(IOUT, "\nHistoric catches:\n");
  137. 	TOTALC = 0;
  138. 	for (I = 0; I <= MAXYR; I++) {
  139. 		_getline(buff, sizeof(buff), IN);
  140. 		sscanf(buff,"%4d%7lf", &IY, &C);
  141. 		if (IY < 0) break; /* EOD */
  142.  
  143. 		fprintf(IOUT, "%4d%6.f\n", IY, C);
  144. 		IY = IY - INITYR;
  145. 		if (IY < 0 || IY > md.IYEAR) {
  146. 			fprintf(stderr, " ERROR: IY is out of range(%d)\n", IY);
  147. 			ret = 1;
  148. 			goto exit;
  149. 		}
  150. 		CATCH[IY] = C;
  151. 		TOTALC = TOTALC + C;
  152. 	}
  153.  
  154. 	if (TOTALC < 0) {
  155. 		fprintf(stderr, " ERROR: No historic catch input\n");
  156. 		ret = 1;
  157. 		goto exit;
  158. 	}
  159.  
  160. /* Read in the non-zero sightings estimates and information matrix */
  161.  
  162. 	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
  163. 	_getline(buff, sizeof(buff), IN); /* Discard empty lines */
  164. 	_getline(buff, sizeof(buff), IN);
  165. 	numlen = 4; len = strlen(buff);
  166. 	md.NS = atoi(&buff[len-numlen]);
  167.  
  168. 	fprintf(IOUT,   "\nAbundance estimates:\n");
  169. 	if (md.NS > MAXEST) {
  170. 		fprintf(stderr, "ERROR: Abundance year out of range\n");
  171. 		ret = 1;
  172. 		goto exit;
  173. 	}
  174. 	_getline(buff, sizeof(buff), IN);
  175. 	for (N = 0; N <= md.NS - 1; N++) {
  176. 		N1 = (N * (N + 1)) / 2;
  177. 		_getline(buff, sizeof(buff), IN);
  178. 		sscanf(buff, "%4d%10lf", &ISYR[N], &SIGHT[N]);
  179. 		fprintf(IOUT,   "%4d %8.f", ISYR[N], SIGHT[N]);
  180. 		len = strlen(buff);
  181. 		for (J = 0; J < N + 1; J++) {
  182. 			if (len > J * 7 + 14) {
  183. 				sscanf(&buff[0] + J * 7 + 14,"%7lf", &FMATRX[N1 + J]);
  184. 			} else {
  185. 				FMATRX[N1 + J] = 0.0;
  186. 			}
  187. 			fprintf(IOUT,   "%7.f", FMATRX[N1 + J]);
  188. 		}
  189. 		fprintf(IOUT,   "\n");
  190. 		ISYR[N] = ISYR[N] - INITYR;
  191. 		if (ISYR[N] < 0 || ISYR[N] >= md.IYEAR) {
  192. 			fprintf(stderr, "ERROR: Sight year out of range\n");
  193. 			ret = 1;
  194. 			goto exit;
  195. 		}
  196. 		if (SIGHT[N] <= 0.0) {
  197. 			fprintf(stderr, "ERROR: Estimate not positive\n");
  198. 			ret = 1;
  199. 			goto exit;
  200. 		}
  201. 	}
  202.  
  203. /* Read in the Poisson multipliers for any zero sightings estimates */
  204.  
  205. 	_getline(buff, sizeof(buff), IN);
  206. 	_getline(buff, sizeof(buff), IN);
  207. 	numlen = 3; len = strlen(buff);
  208. 	md.NZ = atoi(&buff[len-numlen]);
  209. 	if (md.NZ > 0) {
  210. 		fprintf(IOUT,   "\nZero abundance estimates:\n");
  211. 	}
  212. 	if (md.NZ > MAXEST) {
  213. 		fprintf(stderr, "ERROR: Zero estimate array too small\n");
  214. 		ret = 1;
  215. 		goto exit;
  216. 	}
  217. 	_getline(buff, sizeof(buff), IN);
  218. 	for (N = 0; N <= md.NZ - 1; N++) {
  219. 		N1 = (N * (N + 1)) / 2;
  220. 		_getline(buff, sizeof(buff), IN);
  221. 		sscanf(buff, "%4d%10lf", &IZYR[N], &ZMULT[N]);
  222. 		fprintf(IOUT,   "%4d%7.f\n", IZYR[N], ZMULT[N]);
  223. 		IZYR[N] = IZYR[N] - INITYR;
  224. 		if (IZYR[N] < 0 || IZYR[N] >= md.IYEAR) {
  225. 			fprintf(stderr, "Sight year out of range\n");
  226. 			ret = 1;
  227. 			goto exit;
  228. 		}
  229.  
  230. 		if (ZMULT[N] <= 0.0) {
  231. 			fprintf(stderr, "ERROR: Multiplier not positive\n");
  232. 			ret = 1;
  233. 			goto exit;
  234. 		}
  235. 	}
  236. 	fprintf(IOUT,"\n");
  237. /* Bound the range for the integration over K */
  238.  
  239. 	md.RKHI = 1.0E7;
  240. 	md.RKLO = 0.0;
  241. 	IS = 0; /* orignal not initialized bug? */
  242.  
  243. /* ****** Run the CLA to obtain the nominal catch limit */
  244. 	RAWCL = CLIMIT(CATCH, SIGHT, FMATRX, ISYR, IZYR, ZMULT, POP, G);
  245. /* Set the catch limits for PCYCLE years. If the catch limit may be */
  246. /* subject to phaseout, call POUT to apply the phaseout rule. */
  247.  
  248. /* First set ILAST = year of the most recent abundance estimate */
  249.  
  250. 	if (md.NS > 0) ILAST = ISYR[md.NS - 1];
  251. 	if (md.NZ > 0) ILAST = (IS > IZYR[md.NZ - 1])? IS: IZYR[md.NZ - 1];
  252.  
  253. 	for (IY = md.IYEAR; IY <= md.IYEAR + mp.PCYCLE - 1; IY++) {
  254. 		if (POUT == 1) {
  255. 			PHOUT(&CL, RAWCL, ILAST, IY);
  256. 		} else {
  257. 			CL = RAWCL;
  258. 		}
  259. 		fprintf(IOUT,   "Year: %4d     Catch limit:%6d\n", IY + INITYR, (int)(CL + 0.5));
  260.  
  261. 	}
  262.  
  263. exit:
  264. 	return ret;
  265. }
  266.  
  267. /*
  268. C **************************************************************************************
  269. C
  270. C	CLC version 6
  271. C
  272. C **************************************************************************************
  273. C
  274. C	31 January 1994
  275. C
  276. C **************************************************************************************
  277. */
  278. /*
  279. 	SUBROUTINE PHOUT (CL, RAWCL, ILAST, IY)
  280. */
  281. void PHOUT(double *CL, double RAWCL, int ILAST, int IY)
  282. {
  283. /*
  284. 	PHASET Number of years without surveys before phaseout invoked
  285. 	PHASEP Phaseout annual reduction proportion
  286. 	Phaseout: Reduce catch limit if there is no survey data in the
  287. 	last PHASET years
  288. */
  289. 	if (IY >= ILAST + mp.PHASET) {
  290. 		*CL = RAWCL * (1.0 - mp.PHASEP * (IY - ILAST - mp.PHASET));
  291. 		if (*CL < 0.0) *CL = 0.0;
  292. 	} else {
  293. 		*CL = RAWCL;
  294. 	}
  295. 	/* RETURN
  296. 	   END */
  297. }
  298.  
  299. /*
  300. 	REAL FUNCTION CLIMIT (CATCH,SIGHT,FMATRX,ISYR,IZYR,ZMULT,POP,G)
  301. 	Run the CLA to obtain to nominal catch limit
  302. */
  303. double CLIMIT(double *CATCH, double *SIGHT, double *FMATRX, int *ISYR, int *IZYR, double *ZMULT,double *POP, double *G)
  304. {
  305. 	double SF, SS0, SS1, SS2, SS3;
  306. 	double Y[MAXSTP + 1], B[MAXSTP + 1], RLGB[MAXSTP + 1];
  307. 	double *PRES = NULL, *QRES = NULL;
  308.  
  309. 	int NB, NR;
  310. 	double BINC, YINC, DK, R, RK, D, DD, P, Q, PTOT;
  311. 	int I, J, N, N2;
  312.  
  313. 	double _CLIMIT = 0.0;
  314. 	if (md.NS <= 0) return _CLIMIT;
  315.  
  316. 	PRES = (double *) malloc(sizeof(double) * (MAXSIM + 1));
  317. 	if (PRES == NULL) {
  318. 		fprintf(stderr, "Can't Allocate memory (PRES)\n");
  319. 		goto err_exit;
  320. 	}
  321. 	memset(PRES, 0, sizeof(double) * (MAXSTP + 1));
  322.  
  323. 	QRES = (double *) malloc(sizeof(double) * (MAXSIM + 1));
  324. 	if (QRES == NULL) {
  325. 		fprintf(stderr, "Can't Allocate memory (QRES)\n");
  326. 		goto err_exit;
  327. 	}
  328. 	memset(QRES, 0, sizeof(double) * (MAXSTP + 1));
  329.  
  330. 	SF = 0.5 / (mp.PSCALE * mp.PSCALE);
  331.  
  332. 	/* Check the sizes of the Y and B arrays are large enough */
  333. 	if (mp.PNBSTP > MAXSTP || mp.PNYSTP > MAXSTP) {
  334. 		fprintf(stderr, "Y &/or b array sizes too small\n");
  335. 		goto err_exit;
  336. 	}
  337. 	/* Set sightings bias step sizes & their log values. BINC = Bias increment */
  338. 	NB = (int)mp.PNBSTP;
  339. 	BINC = (mp.PBMAX - mp.PBMIN) / NB;
  340. 	for (I = 0; I <= NB - 1; I++) {
  341.         *(B + I) = mp.PBMIN + (I + 0.5) * BINC;
  342.         *(RLGB + I) = -log(*(B + I));
  343. 	}
  344.  
  345. 	/* Set productivity parameter step sizes (midpoints) */
  346. 	NR = (int)mp.PNYSTP;
  347. 	YINC = mp.PYMAX / NR;
  348. 	for (I = 0; I <= NR - 1; I++) {
  349. 		*(Y + I) = (I + 0.5) * YINC;
  350. 	}
  351.  
  352. 	PTOT = 0;
  353. 	N = 0;
  354. 	for (I = 0; I <= NR - 1; I++) {
  355. 		/* Set R from the productivity parameter, Y */
  356. 		R = 1.4184 * *(Y + I);
  357.  
  358. 		/* Step size for K */
  359.         DK = mp.PKSTEP;
  360.         D = 1.0;
  361.         RK = md.RKHI;
  362.  
  363. 		/* Use function STKSIM to set up the Nth population trajectory */
  364. 		/* i.e. set up the pop array */
  365.  
  366. 		while (1) {
  367. 			if ((RK <= md.RKLO) || (STKSIM(RK, R, POP, CATCH) <= 0.0)) break;
  368. 			if (N >= MAXSIM) {
  369. 				fprintf(stderr, "ERROR: TOO MANY SIMULATIONS");
  370. 				goto err_exit;
  371. 			}
  372. 			/* How much depletion covered? */
  373. 			DD = D - POP[md.IYEAR] / RK;
  374. 			D = POP[md.IYEAR] / RK;
  375. 			P = 0.0;
  376.  
  377. 			if (DD > 0.0) {
  378. 				/* Compute the internal catch limit corresponding to D and Y(I) */
  379. 				QRES[N] = CONTRL(D, Y[I], mp.PLEVEL, mp.PSLOPE) * POP[md.IYEAR];
  380.  
  381. 				/* Calculate deviance */
  382.             	DEVIAN(&SS0, &SS1, &SS2, &SS3, SIGHT, FMATRX, ISYR,IZYR, ZMULT, POP, G);
  383.  
  384. 				/* Scale likelihood and integrate over values for the bias parameter */
  385. 				for (J = 0; J <= NB - 1; J++) {
  386. 					P = P + exp(-SF * (SS0 + *(RLGB + J) * (SS1 + *(RLGB + J) * SS2) + SS3 * *(B + J)));
  387. 				}
  388.  
  389. 				/* Calculate the weight for this point (& total area under likelihood) */
  390. 				PRES[N] = P * DD;
  391. 				PTOT = PTOT + PRES[N];
  392.  
  393. 				/* Update counter */
  394. 				N = N + 1;
  395.  
  396. 				/* Find the next K */
  397. 				DK = DK * mp.PDSTEP / DD;
  398. 				if (DK > mp.PKSTEP) DK = mp.PKSTEP;
  399.  
  400. 			} else {
  401. 				DK = mp.PKSTEP;
  402. 			}
  403. 			/* Set the new value of K */
  404. 			RK = RK / (1.0 + DK);
  405. 		}
  406. 	}
  407. 	if (PTOT < 0.0) {
  408. 		fprintf(stderr,"ERROR: PROB INTEGRATES TO ZERO");
  409. 		goto err_exit;
  410. 	}
  411.  
  412. 	/* Sort the QRES and PRES arrays in ascending order of QRES */
  413. 	N2 = N;
  414. 	SORT(QRES, PRES, N2);
  415.  
  416. 	/* Normalize the relative likelihoods */
  417. 	for (I = 0; I <= N - 1; I++) {
  418. 		PRES[I] = PRES[I] / PTOT;
  419. 	}
  420.  
  421. 	/* Extract the desired probability level: the nominal catch limit (NCL) */
  422. 	/* is the lower PROB% of the distribution. */
  423. 	/* First calculate PRES(I), the probability that the NCL is between */
  424. 	/* QRES(I) & QRES(I + 1). */
  425.  
  426. 	P = 0;
  427. 	for (I = 0; I <= N - 1; I++) {
  428. 		P = P + PRES[I];
  429. 		if (P > mp.PPROB) break;
  430. 	}
  431.  
  432. 	/* Interpolate to set the nominal catch limit */
  433. 	if  (I >= N - 1) {
  434.         Q = QRES[N - 1];
  435. 	} else {
  436.         Q = (QRES[I + 1] * (mp.PPROB - P + PRES[I]) + QRES[I] * (P - mp.PPROB)) / PRES[I];
  437. 	}
  438.  
  439. 	_CLIMIT = Q;
  440.  
  441. 	if (QRES != NULL) free(QRES);
  442. 	if (PRES != NULL) free(PRES);
  443. 	return _CLIMIT;
  444.  
  445. err_exit:
  446. 	if (QRES != NULL) free(QRES);
  447. 	if (PRES != NULL) free(PRES);
  448. 	exit(1);
  449. 	return 0.0; /* not reached */
  450. }
  451.  
  452. /*
  453. 	FUNCTION CONTRL (D, Y, PLEVEL, PSLOPE)
  454. 	Catch control law
  455. */
  456. double CONTRL(double D, double Y, double PLEVEL, double PSLOPE)
  457. {
  458. 	double ret;
  459. 	if (D < PLEVEL) {
  460. 		ret = 0.0;
  461. 	} else if (D < 1.0) {
  462. 		ret = PSLOPE * Y * (D - PLEVEL);
  463. 	} else {
  464. 		ret = PSLOPE * Y * (1.0 - PLEVEL);
  465. 	}
  466. 	return ret;
  467. 	/* END */
  468. }
  469.  
  470. /*
  471. 	SUBROUTINE DEVIAN (SS0, SS1, SS2, SS3, SIGHT, FMATRX, ISYR,
  472. 	Calculate deviance (-2 log likelihood) in terms of coefficients for 
  473. 	the bias and log bias parameters
  474. */
  475. void DEVIAN(double *SS0, double *SS1, double *SS2, double *SS3, double *SIGHT,
  476.             double *FMATRX, int *ISYR, int *IZYR, double *ZMULT, double *POP, double *G)
  477. {
  478. 	int N, J, K;
  479.  
  480. 	*SS0 = 0.0;
  481. 	*SS1 = 0.0;
  482. 	*SS2 = 0.0;
  483. 	*SS3 = 0.0;
  484.  
  485. 	for (N = 0; N <= md.NS - 1; N++) { 
  486.  
  487. 		*(G + N) = log(*(SIGHT + N) / *(POP+*(ISYR + N)));
  488. 		K = N * (N + 1) / 2;
  489.  
  490. 		for (J = 0; J <= N - 1; J++) {
  491. 			/*1st add non diagonal contributions (which are doubled up)*/
  492. 			*SS0 = *SS0 + 2.0 * *(G + J) * *(G + N) * *(FMATRX + K);
  493. 			*SS1 = *SS1 + 2.0 * (*(G + J) + *(G + N)) * *(FMATRX + K);
  494. 			*SS2 = *SS2 + *(FMATRX + K) + *(FMATRX + K);
  495. 			K = K + 1;
  496. 		}
  497.  
  498. 		/* Now add diagnoal contribution */
  499.         *SS0 = *SS0 + *(G + N) * *(G + N) * *(FMATRX + K);
  500.         *SS1 = *SS1 + 2.0 * *(G + N) * *(FMATRX + K);
  501.         *SS2 = *SS2 + *(FMATRX + K);
  502.  
  503. 	}
  504.  
  505. 	/* Now do the zero estimates */
  506. 	for (N =0; N <= md.NZ - 1; N++) {
  507. 		*SS3 = *SS3 + 2.0 * *(POP + *(IZYR + N)) / *(ZMULT + N);
  508.  	}
  509.  
  510. 	return;
  511. }
  512. /*
  513. 	FUNCTION STKSIM (RK, R, POP, CATCH)
  514. 	Calculate the stock trajectory with parameter RK and R: return
  515. 	the current stock size
  516. 	RK = notional carrying capacity; R = productivity parameter * 1.4184 
  517. */
  518. double STKSIM(double RK, double R, double *POP, double *CATCH)
  519. {
  520. 	int I;
  521. 	double D;
  522. 	*(POP+md.ISTART) = RK;
  523. 	for (I = md.ISTART + 1; I <= md.IYEAR; I++) {
  524. 		D = *(POP + I - 1) / RK;
  525. 		*(POP + I) = *(POP + I - 1) * (1.0 + R * (1.0 - D * D)) - *(CATCH + I -1);
  526. 		if (*(POP + I) <= 0) {
  527. 			return 0.0;
  528. 		}
  529. 	}
  530.  
  531. 	return *(POP + md.IYEAR);
  532. }
  533. /* 
  534. 	SUBROUTINE SORT (ARRAY, ARRAY2, N)
  535. 	SORT sorts a pair of arrays in ascending order of the first array
  536.    (using the heapsort algorithm)
  537. */
  538. void SORT(double *ARRAY, double *ARRAY2, int N)
  539. {
  540. 	double /*ARRAY(0:*), ARRAY2(0:*),*/ TEMP, TEMP2;
  541. 	int  /*N,*/ K, IR, I, J;
  542. 	if (N < 2) return;
  543. 	K = N / 2;
  544. 	IR = N - 1;
  545. 	while (1) { /*10*/
  546. 		if (K != 0) {
  547. 			K = K - 1;
  548. 			TEMP = *(ARRAY + K);
  549. 			TEMP2 = *(ARRAY2 + K);
  550. 		} else {
  551. 			TEMP = *(ARRAY + IR);
  552. 			TEMP2 = *(ARRAY2 + IR);
  553. 			*(ARRAY + IR) = *ARRAY;
  554. 			*(ARRAY2 + IR) = *ARRAY2;
  555. 			IR = IR - 1;
  556. 			if (IR == 0) {
  557. 				*ARRAY = TEMP;
  558. 				*ARRAY2 = TEMP2;
  559. 				return;
  560. 			}
  561. 		}
  562. 		I = K;
  563. 		J = K + K + 1;
  564. 		while (1) {/*20*/
  565. 			if (J > IR) break;
  566. 			if ((J < IR) && (*(ARRAY + J) < *(ARRAY+ J + 1))) J = J + 1;
  567. 			if (TEMP < *(ARRAY + J)) {
  568. 				*(ARRAY + I) = *(ARRAY + J);
  569. 				*(ARRAY2 + I) = *(ARRAY2 + J);
  570. 				I = J;
  571. 				J = J + I + 1;
  572. 			} else {
  573. 				J = IR + 1;
  574. 			}
  575. 		} /*GOTO 20*/
  576. 		*(ARRAY + I) = TEMP;
  577. 		*(ARRAY2 + I) = TEMP2;
  578. 	}/*GOTO 10*/
  579. 	/*END*/
  580. }
  581. int _getline(char *buff, size_t size, FILE *fp){
  582. 	char *ret;
  583. 	ret = fgets(buff, size, fp);
  584. 	if (ret != NULL) {
  585. 		if (buff[strlen(buff) - 1] == '\n')
  586. 			buff[strlen(buff) - 1] = 0;
  587. 		if (buff[strlen(buff) - 1] == '\r')
  588. 			buff[strlen(buff) - 1] = 0;
  589. 	}
  590. 	return (ret != NULL);
  591. }
  592.  
Success #stdin #stdout 0.18s 9432KB
comments (?)
stdin
copy 
TEST STOCK (FOR A SMALL AREA)

First year of input data     1988
Year of next catch limit     1994
Test for phaseout               1

Catch data:
(I4,F8.0)
1988    20.
1989    20.
1990    20.
1991    20.
1992    20.
1993    20.
-1

abundance data:
No. of sightings estimate      1
(I4,F10.0,10F8.0)
1988   27150.0   101.

No. of zero estimate           0
(I4,F10.0)
stdout
copy 
TEST STOCK (FOR A SMALL AREA) by C

First year of input data      1988
Year of next catch limit      1994
Apply phaseout if necessary    YES

Historic catches:
1988    20
1989    20
1990    20
1991    20
1992    20
1993    20

Abundance estimates:
1988    27150    101

Year: 1994     Catch limit:   493
Year: 1995     Catch limit:   493
Year: 1996     Catch limit:   493
Year: 1997     Catch limit:   395
Year: 1998     Catch limit:   296
https://ideone.com/uoOEqd
language:
C (gcc 8.3)
created:
5 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!