/*
Copyright (C) 1995 Evgenii Rudnyi, http://Evgenii.Rudnyi.Ru

Software for paper 
E.B. Rudnyi. Statistical model of systematic errors: linear error model. Chemometrics and Intelligent Laboratory Systems. 1996, V. 34, N 1, p. 41-54. 

This software is a copyrighted work licensed under the terms, described
in the file "FREE_LICENSE". 
*/

/* hscor.f -- translated by f2c (version 19950211).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include <iostream.h>
#include <iomanip.h>

#ifdef __cplusplus
extern "C" {
#endif
#include "g2c.h"

/* Common Block Declarations */

struct {
    doublereal hmelt, tmelt, cs[5], cl[5], cm[5], cd[5];
} hsc_
= {26320., 1044.,
//KCl(s) 298.15 - 1044
  0.440080000000E+02, -0.200000000000E-05,  0.331650000000E+03,
 -0.359910000000E+04,  0.228500000004E+05,
//KCl(l) 1044 - 3000
  0.720000000000E+02,  0.000000000000E+00,  0.000000000000E+00,
  0.000000000000E+00,  0.000000000000E+00,
//KCl(g) 298.15 - 1500
  0.371655807490E+02, -0.909746275280E-03,  0.131289806366E+02,
 -0.374903526306E+02,  0.151497642516E+03,
//K2Cl2(g) 298.15 - 1500
  0.829382476810E+02, -0.208506127820E-02,  0.445766353600E+01,
 -0.350137882230E+02,  0.942856521600E+02};

#define hsc_1 hsc_

/*<       SUBROUTINE HSCOR (TF, DH1, DS1, DH2, DS2) >*/
/* Subroutine */ int hscor_(doublereal *tf, doublereal *dh1, doublereal *ds1, 
	doublereal *dh2, doublereal *ds2)
{
    /* System generated locals */
    doublereal d__1, d__2, d__3, d__4, d__5, d__6, d__7, d__8;

    /* Builtin functions */
    double log(doublereal);

    /* Local variables */
    doublereal t, smelt, t1;

/*<       IMPLICIT DOUBLE PRECISION (A-H, O-Z) >*/
/*<       DIMENSION CS(5), CL(5), CM(5), CD(5) >*/
/*<       COMMON /HSC/ HMELT, TMELT, CS, CL, CM, CD >*/
/*<       SMELT=HMELT/TMELT >*/
    smelt = hsc_1.hmelt / hsc_1.tmelt;
/*<       T1=TMELT*1E-4 >*/
    t1 = hsc_1.tmelt * 1e-4;
/*<       T=TF*1E-4 >*/
    t = *tf * 1e-4;
/*<       IF (T.LT.T1) GO TO 5 >*/
    if (t < t1) {
	goto L5;
    }
/*<       IF (T.GT.T1) GO TO 10 >*/
    if (t > t1) {
	goto L10;
    }
/*<       DH1=0. >*/
    *dh1 = 0.;
/*<       DS1=0. >*/
    *ds1 = 0.;
/*<       DH2=0. >*/
    *dh2 = 0.;
/*<       DS2=0. >*/
    *ds2 = 0.;
/*<       RETURN >*/
    return 0;
/*< 5  >*/
L5:
/* Computing 2nd power */
    d__1 = t;
/* Computing 2nd power */
    d__2 = t1;
/* Computing 3rd power */
    d__3 = t, d__4 = d__3;
/* Computing 3rd power */
    d__5 = t1, d__6 = d__5;
/* Computing 4th power */
    d__7 = t, d__7 *= d__7;
/* Computing 4th power */
    d__8 = t1, d__8 *= d__8;
    *dh1 = ((hsc_1.cm[0] - hsc_1.cs[0]) * (t - t1) - (hsc_1.cm[1] - hsc_1.cs[
            1]) * (1. / t - 1. / t1) + (hsc_1.cm[2] - hsc_1.cs[
            2]) / 2. * (d__1 * d__1 - d__2 * d__2) + (hsc_1.cm[3] -
            hsc_1.cs[3]) / 3. * (d__4 * (d__3 * d__3) - d__6 * (d__5 *
            d__5)) + (hsc_1.cm[4] - hsc_1.cs[4]) * .25 * (d__7 * d__7
            - d__8 * d__8)) * 1e4;
/*<    >*/
/* Computing 2nd power */
    d__1 = t;
/* Computing 2nd power */
    d__2 = t1;
/* Computing 2nd power */
    d__3 = t;
/* Computing 2nd power */
    d__4 = t1;
/* Computing 3rd power */
    d__5 = t, d__6 = d__5;
/* Computing 3rd power */
    d__7 = t1, d__8 = d__7;
    *ds1 = (hsc_1.cm[0] - hsc_1.cs[0]) * log(t / t1) - (hsc_1.cm[1] - 
            hsc_1.cs[1]) * .5 * (1. / (d__1 * d__1) - 1.
	    / (d__2 * d__2)) + (hsc_1.cm[2] - hsc_1.cs[2]) * (t - t1) + (
            hsc_1.cm[3] - hsc_1.cs[3]) / 2. * (d__3 * d__3 - d__4 *
            d__4) + (hsc_1.cm[4] - hsc_1.cs[4]) / 3. * (d__6 * (d__5 *
	    d__5) - d__8 * (d__7 * d__7));
/*<    >*/
/* Computing 2nd power */
    d__1 = t;
/* Computing 2nd power */
    d__2 = t1;
/* Computing 3rd power */
    d__3 = t, d__4 = d__3;
/* Computing 3rd power */
    d__5 = t1, d__6 = d__5;
/* Computing 4th power */
    d__7 = t, d__7 *= d__7;
/* Computing 4th power */
    d__8 = t1, d__8 *= d__8;
    *dh2 = ((hsc_1.cd[0] - hsc_1.cs[0] * 2.) * (t - t1) - (hsc_1.cd[1]
            - hsc_1.cs[1] * 2.) * (1. / t - 1. / t1) + (
            hsc_1.cd[2] - hsc_1.cs[2] * 2.) / 2. * (d__1 * d__1
            - d__2 * d__2) + (hsc_1.cd[3] - hsc_1.cs[3] * 2.) /
	    3. * (d__4 * (d__3 * d__3) - d__6 * (d__5 * d__5)) + (hsc_1.cd[4] 
            - hsc_1.cs[4] * 2.) * .25 * (d__7 * d__7 - d__8 *
            d__8)) * 1e4;
/*<    >*/
/* Computing 2nd power */
    d__1 = t;
/* Computing 2nd power */
    d__2 = t1;
/* Computing 2nd power */
    d__3 = t;
/* Computing 2nd power */
    d__4 = t1;
/* Computing 3rd power */
    d__5 = t, d__6 = d__5;
/* Computing 3rd power */
    d__7 = t1, d__8 = d__7;
    *ds2 = (hsc_1.cd[0] - hsc_1.cs[0] * 2.) * log(t / t1) - (hsc_1.cd[
            1] - hsc_1.cs[1] * 2.) * .5 * (1. / (d__1 *
            d__1) - 1. / (d__2 * d__2)) + (hsc_1.cd[2] - hsc_1.cs[2] *
            2.) * (t - t1) + (hsc_1.cd[3] - hsc_1.cs[3] * 2.) /
            2. * (d__3 * d__3 - d__4 * d__4) + (hsc_1.cd[4] - hsc_1.cs[
            4] * 2.) / 3. * (d__6 * (d__5 * d__5) - d__8 * (
	    d__7 * d__7));
/*<       DH2 = DH2 - DH1 >*/
    *dh2 -= *dh1;
/*<       DS2 = DS2 - DS1 >*/
    *ds2 -= *ds1;
/*<       RETURN >*/
    return 0;
/*< 10 >*/
L10:
/* Computing 2nd power */
    d__1 = t;
/* Computing 2nd power */
    d__2 = t1;
/* Computing 3rd power */
    d__3 = t, d__4 = d__3;
/* Computing 3rd power */
    d__5 = t1, d__6 = d__5;
/* Computing 4th power */
    d__7 = t, d__7 *= d__7;
/* Computing 4th power */
    d__8 = t1, d__8 *= d__8;
    *dh1 = ((hsc_1.cm[0] - hsc_1.cl[0]) * (t - t1) - (hsc_1.cm[1] - hsc_1.cl[
            1]) * (1. / t - 1. / t1) + (hsc_1.cm[2] - hsc_1.cl[
            2]) / 2. * (d__1 * d__1 - d__2 * d__2) + (hsc_1.cm[3] -
            hsc_1.cl[3]) / 3. * (d__4 * (d__3 * d__3) - d__6 * (d__5 *
            d__5)) + (hsc_1.cm[4] - hsc_1.cl[4]) * .25 * (d__7 * d__7
            - d__8 * d__8)) * 1e4 - hsc_1.hmelt;
/*<    >*/
/* Computing 2nd power */
    d__1 = t;
/* Computing 2nd power */
    d__2 = t1;
/* Computing 2nd power */
    d__3 = t;
/* Computing 2nd power */
    d__4 = t1;
/* Computing 3rd power */
    d__5 = t, d__6 = d__5;
/* Computing 3rd power */
    d__7 = t1, d__8 = d__7;
    *ds1 = (hsc_1.cm[0] - hsc_1.cl[0]) * log(t / t1) - (hsc_1.cm[1] - 
            hsc_1.cl[1]) * .5 * (1. / (d__1 * d__1) - 1.
	    / (d__2 * d__2)) + (hsc_1.cm[2] - hsc_1.cl[2]) * (t - t1) + (
            hsc_1.cm[3] - hsc_1.cl[3]) / 2. * (d__3 * d__3 - d__4 *
            d__4) + (hsc_1.cm[4] - hsc_1.cl[4]) / 3. * (d__6 * (d__5 *
	    d__5) - d__8 * (d__7 * d__7)) - smelt;
/*<    >*/
/* Computing 2nd power */
    d__1 = t;
/* Computing 2nd power */
    d__2 = t1;
/* Computing 3rd power */
    d__3 = t, d__4 = d__3;
/* Computing 3rd power */
    d__5 = t1, d__6 = d__5;
/* Computing 4th power */
    d__7 = t, d__7 *= d__7;
/* Computing 4th power */
    d__8 = t1, d__8 *= d__8;
    *dh2 = ((hsc_1.cd[0] - hsc_1.cl[0] * 2.) * (t - t1) - (hsc_1.cd[1]
            - hsc_1.cl[1] * 2.) * (1. / t - 1. / t1) + (
            hsc_1.cd[2] - hsc_1.cl[2] * 2.) / 2. * (d__1 * d__1
            - d__2 * d__2) + (hsc_1.cd[3] - hsc_1.cl[3] * 2.) /
	    3. * (d__4 * (d__3 * d__3) - d__6 * (d__5 * d__5)) + (hsc_1.cd[4] 
            - hsc_1.cl[4] * 2.) * .25 * (d__7 * d__7 - d__8 *
            d__8)) * 1e4 - hsc_1.hmelt * 2.;
/*<    >*/
/* Computing 2nd power */
    d__1 = t;
/* Computing 2nd power */
    d__2 = t1;
/* Computing 2nd power */
    d__3 = t;
/* Computing 2nd power */
    d__4 = t1;
/* Computing 3rd power */
    d__5 = t, d__6 = d__5;
/* Computing 3rd power */
    d__7 = t1, d__8 = d__7;
    *ds2 = (hsc_1.cd[0] - hsc_1.cl[0] * 2.) * log(t / t1) - (hsc_1.cd[
            1] - hsc_1.cl[1] * 2.) * .5 * (1. / (d__1 *
            d__1) - 1. / (d__2 * d__2)) + (hsc_1.cd[2] - hsc_1.cl[2] *
            2.) * (t - t1) + (hsc_1.cd[3] - hsc_1.cl[3] * 2.) /
            2. * (d__3 * d__3 - d__4 * d__4) + (hsc_1.cd[4] - hsc_1.cl[
            4] * 2.) / 3. * (d__6 * (d__5 * d__5) - d__8 * (
            d__7 * d__7)) - smelt * 2.;
/*<       DH2 = DH2 - DH1 >*/
    *dh2 -= *dh1;
/*<       DS2 = DS2 - DS1 >*/
    *ds2 -= *ds1;
/*<       RETURN >*/
    return 0;
/*<       END >*/
} /* hscor_ */

#ifdef __cplusplus
	}
#endif

/*
int main()
{
  double dh1, ds1, dh2, ds2;
  for (double T = 800.; T <= 1601.; T +=100)
  {
    hscor_(&T, &dh1, &ds1, &dh2, &ds2);
    cout << setw(10) << T
         << setw(10) << dh1
         << setw(10) << ds1
         << setw(10) << dh2
         << setw(10) << ds2
         << endl;
  }
}
*/