SCons can now build the contents of the 'ext' directory
This commit is contained in:
Ray Speth
parent
5345040ec0
commit
2c4623e336
30 changed files with 134 additions and 1423 deletions
70
SConstruct
70
SConstruct
|
|
@ -50,28 +50,32 @@ opts.AddVariables(
|
|||
BoolVariable('lapack_ftn_trailing_underscore', '', True),
|
||||
BoolVariable('lapack_ftn_string_len_at_end', '', True),
|
||||
('bitcompile', '', ''), # '32' or '64'
|
||||
('cxx', '', 'g++'),
|
||||
('cc', '', 'gcc'),
|
||||
('cxxflags', '', '-O3 -Wall'),
|
||||
('lcxx_end_libs', '-lm'),
|
||||
('pic', '', ''),
|
||||
('shared', '', '-dynamic'),
|
||||
('CXX', '', env['CXX']),
|
||||
('CC', '', env['CC']),
|
||||
('CXXFLAGS', '', '-O3 -Wall'),
|
||||
BoolVariable('build_thread_safe', '', False),
|
||||
BoolVariable('build_with_f2c', '', True),
|
||||
('f77', '', 'g77'),
|
||||
('fflags', '', '-O3'),
|
||||
('lfort_flags', '', '-L/usr/local/lib'),
|
||||
('archive', '', 'ar ruv'),
|
||||
('ranlib', '', 'ranlib'),
|
||||
('F77', '', env['F77']),
|
||||
('F77FLAGS', '', '-O3'),
|
||||
('F90', '', env['F90']),
|
||||
('F90FLAGS', '', '-O3'),
|
||||
('install_bin', '', 'config/install-sh'),
|
||||
('graphvisdir', '' ,''),
|
||||
('cxx_ext', '', 'cpp'),
|
||||
('f77_ext', '', 'f'),
|
||||
('f90_ext', '', 'f90'),
|
||||
('exe_ext', '', ''),
|
||||
('ct_shared_lib', '', 'clib'),
|
||||
('rpfont', '', 'Helvetica'),
|
||||
('cantera_version', '', '1.8.x')
|
||||
# These variables shouldn't be necessary any more...
|
||||
# ('cxx_ext', '', 'cpp'),
|
||||
# ('f77_ext', '', 'f'),
|
||||
# ('f90_ext', '', 'f90'),
|
||||
# ('exe_ext', '', ''),
|
||||
# ('lcxx_end_libs', '-lm'),
|
||||
# ('pic', '', '-fPIC'),
|
||||
# ('shared', '', '-dynamic'),
|
||||
# ('lfort_flags', '', '-L/usr/local/lib'),
|
||||
# ('AR', '', env['AR']),
|
||||
# ('ARFLAGS', '', env['ARFLAGS']), # ('archive', '', 'ar ruv'),
|
||||
# ('ranlib', '', 'ranlib'),
|
||||
)
|
||||
|
||||
opts.Update(env)
|
||||
|
|
@ -88,13 +92,12 @@ if env['python_package'] in ('full', 'default'):
|
|||
None, PathVariable.PathAccept),
|
||||
PathVariable('cantera_python_home', 'where to install the python package',
|
||||
None, PathVariable.PathAccept),
|
||||
|
||||
)
|
||||
|
||||
# Options that apply only if building the Matlab interface
|
||||
if env['matlab_toolbox'] != 'n':
|
||||
opts.AddVariables(
|
||||
PathVariable('matlab_cmd', 'Path to the matlab executable',
|
||||
PathVariable('matlab_cmd', 'Path to the matlab executable',
|
||||
'default', PathVariable.PathAccept)
|
||||
)
|
||||
|
||||
|
|
@ -125,9 +128,20 @@ opts.Save('cantera.conf', env)
|
|||
# ********************************************
|
||||
env['OS'] = platform.system()
|
||||
|
||||
#def ArithCheck(context):
|
||||
# context.Message('Trying to generate arith.h\n')
|
||||
# exitStatus = context.TryLink(file('ext/f2c_libs/arithchk.c').read(), '.c')
|
||||
# print exitStatus
|
||||
# exitStatus, output = context.TryRun(file('ext/f2c_libs/arithchk.c').read(), '.c')
|
||||
# print exitStatus, output
|
||||
# context.Result(output)
|
||||
# return exitStatus
|
||||
|
||||
conf = Configure(env)
|
||||
#conf = Configure(env, custom_tests = {'ArithCheck':ArithCheck})
|
||||
|
||||
env['HAS_SSTREAM'] = conf.CheckCXXHeader('sstream', '<>')
|
||||
#conf.ArithCheck()
|
||||
|
||||
env = conf.Finish()
|
||||
|
||||
|
|
@ -151,8 +165,8 @@ cdefine('PURIFY_MODE', 'purify')
|
|||
|
||||
# Need to test all of these to see what platform.system() returns
|
||||
configh['SOLARIS'] = 1 if env['OS'] == 'Solaris' else None
|
||||
configh['DARWIN'] = 1 if env['OS'] == 'Darwin' else None
|
||||
configh['CYGWIN'] = 1 if env['OS'] == 'Cygwin' else None
|
||||
configh['DARWIN'] = 1 if env['OS'] == 'Darwin' else None
|
||||
configh['CYGWIN'] = 1 if env['OS'] == 'Cygwin' else None
|
||||
configh['WINMSVC'] = 1 if env['OS'] == 'Windows' else None
|
||||
cdefine('NEEDS_GENERIC_TEMPL_STATIC_DECL', 'OS', 'Solaris')
|
||||
|
||||
|
|
@ -191,8 +205,24 @@ cdefine('THREAD_SAFE_CANTERA', 'build_thread_safe')
|
|||
cdefine('HAS_SSTREAM', 'HAS_SSTREAM')
|
||||
configh['CANTERA_DATA'] = quoted(os.path.join(env['prefix'], 'data'))
|
||||
|
||||
env.AlwaysBuild(env.Command('config.h', 'config.h.in.scons', ConfigBuilder(configh)))
|
||||
config_h = env.Command('config.h', 'config.h.in.scons', ConfigBuilder(configh))
|
||||
#env.AlwaysBuild(config_h)
|
||||
|
||||
# **********************************************
|
||||
# *** Set additional configuration variables ***
|
||||
# **********************************************
|
||||
if env['blas_lapack_libs'] == '':
|
||||
# External BLAS/LAPACK were not given, so we need to compile them
|
||||
env['BUILD_BLAS_LAPACK'] = True
|
||||
env['blas_lapack_libs'] = '-lctlapack -lctblas'
|
||||
|
||||
# *********************
|
||||
# *** Build Cantera ***
|
||||
# *********************
|
||||
build = 'build'
|
||||
env.SConsignFile()
|
||||
env.Append(CPPPATH=os.getcwd())
|
||||
Export('env', 'build', 'config_h')
|
||||
|
||||
VariantDir('build/ext', 'ext', duplicate=0)
|
||||
SConscript('build/ext/SConscript')
|
||||
|
|
|
|||
|
|
@ -1,3 +1,8 @@
|
|||
from glob import glob
|
||||
import os
|
||||
import shutil
|
||||
from os.path import join as pjoin
|
||||
|
||||
class DefineDict(object):
|
||||
def __init__(self, data):
|
||||
self.data = data
|
||||
|
|
@ -12,6 +17,7 @@ class DefineDict(object):
|
|||
else:
|
||||
return '#define %s %s' % (key, self.data[key])
|
||||
|
||||
|
||||
class ConfigBuilder(object):
|
||||
def __init__(self, defines):
|
||||
self.defines = DefineDict(defines)
|
||||
|
|
@ -34,5 +40,44 @@ class ConfigBuilder(object):
|
|||
for key in sorted(self.defines.undefined):
|
||||
print " %-35s %s" % (key, '*undefined*')
|
||||
|
||||
|
||||
class CopyNoPrefix(object):
|
||||
"""
|
||||
Copy a file, ignoring leading directories that are part
|
||||
of 'prefix' (e.g. the variant directory)
|
||||
"""
|
||||
def __init__(self, prefix):
|
||||
self.prefix = prefix
|
||||
|
||||
def __call__(self, source, target, env):
|
||||
sourcepath = psplit(str(source[0]))
|
||||
targetpath = psplit(str(target[0]))
|
||||
|
||||
depth = 0
|
||||
for a,b in zip(targetpath, psplit(self.prefix)):
|
||||
if a == b:
|
||||
depth += 1
|
||||
else:
|
||||
break
|
||||
print str(source[0]), pjoin(*targetpath[depth:])
|
||||
shutil.copyfile(str(source[0]), pjoin(*targetpath[depth:]))
|
||||
|
||||
|
||||
def quoted(s):
|
||||
return '"%s"' % s
|
||||
|
||||
|
||||
def mglob(env, subdir, *args):
|
||||
""" each arg in args is assumed to be file extension """
|
||||
return sum((env.Glob('%s/*.%s' % (subdir, ext)) for ext in args), [])
|
||||
|
||||
|
||||
def psplit(s):
|
||||
head, tail = os.path.split(s)
|
||||
path = [tail]
|
||||
while head:
|
||||
head, tail = os.path.split(head)
|
||||
path.append(tail)
|
||||
|
||||
path.reverse()
|
||||
return path
|
||||
|
|
|
|||
39
ext/SConscript
Normal file
39
ext/SConscript
Normal file
|
|
@ -0,0 +1,39 @@
|
|||
from buildutils import *
|
||||
|
||||
Import('env', 'build', 'config_h')
|
||||
localenv = env.Clone()
|
||||
|
||||
# (subdir, library name, (file extensions))
|
||||
libs = [('tpx','tpx',('cpp')),
|
||||
('math', 'ctmath', ('cpp','c','f'))]
|
||||
|
||||
|
||||
if env['build_with_f2c']:
|
||||
# Create arith.h using the arithchk program
|
||||
arithenv = env.Clone()
|
||||
arithenv.Append(LINKFLAGS='-lm', CPPFLAGS='-DNO_FPINIT') # TODO: make link flag more general
|
||||
arithenv.Program('f2c_libs/arithchk/arithchk', source='f2c_libs/arithchk/arithchk.c')
|
||||
arithenv.Command('arith.h', 'f2c_libs/arithchk/arithchk', '$SOURCE > $TARGET')
|
||||
arithenv.Command('f2c_libs/arith.h', 'arith.h', CopyNoPrefix(build))
|
||||
|
||||
# Possibly system-depenent headers
|
||||
localenv.Command('f2c_libs/signal1.h', 'f2c_libs/signal1.h0', CopyNoPrefix(build))
|
||||
localenv.Command('f2c_libs/sysdep1.h', 'f2c_libs/sysdep1.h0', CopyNoPrefix(build))
|
||||
|
||||
libs.append(('f2c_libs', 'f2c', 'c'))
|
||||
|
||||
if env['BUILD_BLAS_LAPACK']:
|
||||
libs.append(('f2c_blas', 'ctblas', ('c')))
|
||||
libs.append(('f2c_lapack', 'ctlapack', ('c')))
|
||||
|
||||
else:
|
||||
if env['BUILD_BLAS_LAPACK']:
|
||||
libs.append(('blas', 'ctblas', ('f')))
|
||||
libs.append(('lapack', 'ctlapack', ('f')))
|
||||
|
||||
if env['use_sundials'] == 'n':
|
||||
libs.append(('cvode', 'cvode', 'c'))
|
||||
|
||||
for subdir, libname, extensions in libs:
|
||||
lib = localenv.Library(pjoin('../lib', libname),
|
||||
source=mglob(localenv, subdir, *extensions))
|
||||
|
|
@ -1,52 +0,0 @@
|
|||
#include "f2c.h"
|
||||
#include "fio.h"
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
static FILE *
|
||||
#ifdef KR_headers
|
||||
unit_chk(Unit, who) integer Unit; char *who;
|
||||
#else
|
||||
unit_chk(integer Unit, char *who)
|
||||
#endif
|
||||
{
|
||||
if (Unit >= MXUNIT || Unit < 0)
|
||||
f__fatal(101, who);
|
||||
return f__units[Unit].ufd;
|
||||
}
|
||||
|
||||
longint
|
||||
#ifdef KR_headers
|
||||
ftell64_(Unit) integer *Unit;
|
||||
#else
|
||||
ftell64_(integer *Unit)
|
||||
#endif
|
||||
{
|
||||
FILE *f;
|
||||
return (f = unit_chk(*Unit, "ftell")) ? FTELL(f) : -1L;
|
||||
}
|
||||
|
||||
int
|
||||
#ifdef KR_headers
|
||||
fseek64_(Unit, offset, whence) integer *Unit, *whence; longint *offset;
|
||||
#else
|
||||
fseek64_(integer *Unit, longint *offset, integer *whence)
|
||||
#endif
|
||||
{
|
||||
FILE *f;
|
||||
int w = (int)*whence;
|
||||
#ifdef SEEK_SET
|
||||
static int wohin[3] = { SEEK_SET, SEEK_CUR, SEEK_END };
|
||||
#endif
|
||||
if (w < 0 || w > 2)
|
||||
w = 0;
|
||||
#ifdef SEEK_SET
|
||||
w = wohin[w];
|
||||
#endif
|
||||
return !(f = unit_chk(*Unit, "fseek"))
|
||||
|| FSEEK(f, (OFF_T)*offset, w) ? 1 : 0;
|
||||
}
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,39 +0,0 @@
|
|||
#include "f2c.h"
|
||||
//#ifdef __cplusplus
|
||||
//extern "C" {
|
||||
//#endif
|
||||
|
||||
//#ifdef KR_headers
|
||||
//longint pow_qq(ap, bp) longint *ap, *bp;
|
||||
//#else
|
||||
longint pow_qq(longint *ap, longint *bp)
|
||||
//#endif
|
||||
{
|
||||
longint pow, x, n;
|
||||
unsigned long long u; /* system-dependent */
|
||||
|
||||
x = *ap;
|
||||
n = *bp;
|
||||
|
||||
if (n <= 0) {
|
||||
if (n == 0 || x == 1)
|
||||
return 1;
|
||||
if (x != -1)
|
||||
return x == 0 ? 1/x : 0;
|
||||
n = -n;
|
||||
}
|
||||
u = n;
|
||||
for(pow = 1; ; )
|
||||
{
|
||||
if(u & 01)
|
||||
pow *= x;
|
||||
if(u >>= 1)
|
||||
x *= x;
|
||||
else
|
||||
break;
|
||||
}
|
||||
return(pow);
|
||||
}
|
||||
//#ifdef __cplusplus
|
||||
//}
|
||||
//#endif
|
||||
|
|
@ -1,72 +0,0 @@
|
|||
#include "f2c.h"
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
#ifndef LONGBITS
|
||||
#define LONGBITS 32
|
||||
#endif
|
||||
|
||||
#ifndef LONG8BITS
|
||||
#define LONG8BITS (2*LONGBITS)
|
||||
#endif
|
||||
|
||||
longint
|
||||
#ifdef KR_headers
|
||||
qbit_bits(a, b, len) longint a; integer b, len;
|
||||
#else
|
||||
qbit_bits(longint a, integer b, integer len)
|
||||
#endif
|
||||
{
|
||||
/* Assume 2's complement arithmetic */
|
||||
|
||||
ulongint x, y;
|
||||
|
||||
x = (ulongint) a;
|
||||
y = (ulongint)-1L;
|
||||
x >>= b;
|
||||
y <<= len;
|
||||
return (longint)(x & ~y);
|
||||
}
|
||||
|
||||
longint
|
||||
#ifdef KR_headers
|
||||
qbit_cshift(a, b, len) longint a; integer b, len;
|
||||
#else
|
||||
qbit_cshift(longint a, integer b, integer len)
|
||||
#endif
|
||||
{
|
||||
ulongint x, y, z;
|
||||
|
||||
x = (ulongint)a;
|
||||
if (len <= 0) {
|
||||
if (len == 0)
|
||||
return 0;
|
||||
goto full_len;
|
||||
}
|
||||
if (len >= LONG8BITS) {
|
||||
full_len:
|
||||
if (b >= 0) {
|
||||
b %= LONG8BITS;
|
||||
return (longint)(x << b | x >> LONG8BITS - b );
|
||||
}
|
||||
b = -b;
|
||||
b %= LONG8BITS;
|
||||
return (longint)(x << LONG8BITS - b | x >> b);
|
||||
}
|
||||
y = z = (unsigned long)-1;
|
||||
y <<= len;
|
||||
z &= ~y;
|
||||
y &= x;
|
||||
x &= z;
|
||||
if (b >= 0) {
|
||||
b %= len;
|
||||
return (longint)(y | z & (x << b | x >> len - b));
|
||||
}
|
||||
b = -b;
|
||||
b %= len;
|
||||
return (longint)(y | z & (x >> b | x << len - b));
|
||||
}
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,17 +0,0 @@
|
|||
#include "f2c.h"
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
longint
|
||||
#ifdef KR_headers
|
||||
qbit_shift(a, b) longint a; integer b;
|
||||
#else
|
||||
qbit_shift(longint a, integer b)
|
||||
#endif
|
||||
{
|
||||
return b >= 0 ? a << b : (longint)((ulongint)a >> -b);
|
||||
}
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,83 +0,0 @@
|
|||
#/bin/sh
|
||||
#
|
||||
# $Source: /cvsroot/cantera/cantera/ext/f2c_recipes/Makefile.in,v $
|
||||
# $Author: hkmoffa $
|
||||
# $Revision: 1.5 $
|
||||
# $Date: 2008/12/30 21:49:42 $
|
||||
#
|
||||
|
||||
.SUFFIXES :
|
||||
.SUFFIXES : .c .d .o
|
||||
|
||||
do_ranlib = @DO_RANLIB@
|
||||
|
||||
PURIFY=@PURIFY@
|
||||
|
||||
# the directory where the Cantera libraries are located
|
||||
CANTERA_LIBDIR=@buildlib@
|
||||
|
||||
# the directory where Cantera include files may be found.
|
||||
CANTERA_INCDIR=@ctroot@/build/include/cantera
|
||||
|
||||
# the C++ compiler
|
||||
CXX = @CXX@
|
||||
|
||||
# the C compiler
|
||||
CC = @CC@
|
||||
|
||||
# C++ compile flags
|
||||
PIC_FLAG=@PIC@
|
||||
CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT) $(PIC_FLAG)
|
||||
|
||||
# Local include files
|
||||
CXX_INCLUDES=-I../f2c_libs
|
||||
|
||||
# How to compile the dependency file
|
||||
.c.d:
|
||||
@CXX_DEPENDS@ $(CXX_FLAGS) $(CXX_INCLUDES) $*.c > $*.d
|
||||
|
||||
# How to compile a C file
|
||||
.c.o:
|
||||
$(PURIFY) @CC@ -c $< @DEFS@ $(CXX_FLAGS) $(CXX_INCLUDES)
|
||||
|
||||
# -----------------------------------------------
|
||||
|
||||
LIB = $(CANTERA_LIBDIR)/librecipes.a
|
||||
|
||||
all: $(LIB)
|
||||
|
||||
# list of object files
|
||||
OBJS = simp1.o simp2.o simp3.o simplx.o splint.o \
|
||||
splie2.o spline.o splin2.o
|
||||
|
||||
# list of source files
|
||||
SRCS = $(OBJS:.o=.c)
|
||||
|
||||
# List of dependency files to be created
|
||||
DEPENDS=$(OBJS:.o=.d)
|
||||
|
||||
# rule to make library
|
||||
$(LIB): $(OBJS)
|
||||
@ARCHIVE@ $(LIB) $(OBJS) > /dev/null
|
||||
ifeq ($(do_ranlib),1)
|
||||
@RANLIB@ $(LIB)
|
||||
endif
|
||||
|
||||
# ------------------------------------------------
|
||||
# Utility Targets
|
||||
|
||||
clean:
|
||||
$(RM) $(OBJS) $(LIB) *.d .depends
|
||||
|
||||
# depends target
|
||||
depends:
|
||||
$(RM) *.d .depends
|
||||
@MAKE@ .depends
|
||||
|
||||
.depends: $(DEPENDS)
|
||||
cat *.d > .depends
|
||||
|
||||
|
||||
ifeq ($(wildcard .depends), .depends)
|
||||
include .depends
|
||||
endif
|
||||
|
|
@ -1,61 +0,0 @@
|
|||
/* simp1.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Subroutine */ int simp1_(doublereal *a, integer *mp, integer *np, integer *
|
||||
mm, integer *ll, integer *nll, integer *iabf, integer *kp, doublereal
|
||||
*bmax)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer a_dim1, a_offset, i__1;
|
||||
doublereal d__1;
|
||||
|
||||
/* Local variables */
|
||||
static integer k;
|
||||
static doublereal test;
|
||||
|
||||
/* Parameter adjustments */
|
||||
--ll;
|
||||
a_dim1 = *mp;
|
||||
a_offset = 1 + a_dim1;
|
||||
a -= a_offset;
|
||||
|
||||
/* Function Body */
|
||||
*kp = ll[1];
|
||||
*bmax = a[*mm + 1 + (*kp + 1) * a_dim1];
|
||||
if (*nll < 2) {
|
||||
return 0;
|
||||
}
|
||||
i__1 = *nll;
|
||||
for (k = 2; k <= i__1; ++k) {
|
||||
if (*iabf == 0) {
|
||||
test = a[*mm + 1 + (ll[k] + 1) * a_dim1] - *bmax;
|
||||
} else {
|
||||
test = (d__1 = a[*mm + 1 + (ll[k] + 1) * a_dim1], abs(d__1)) -
|
||||
abs(*bmax);
|
||||
}
|
||||
if (test > (float)0.) {
|
||||
*bmax = a[*mm + 1 + (ll[k] + 1) * a_dim1];
|
||||
*kp = ll[k];
|
||||
}
|
||||
/* L11: */
|
||||
}
|
||||
return 0;
|
||||
} /* simp1_ */
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,89 +0,0 @@
|
|||
/* simp2.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Subroutine */ int simp2_(doublereal *a, integer *m, integer *n, integer *
|
||||
mp, integer *np, integer *l2, integer *nl2, integer *ip, integer *kp,
|
||||
doublereal *q1)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer a_dim1, a_offset, i__1, i__2;
|
||||
|
||||
/* Local variables */
|
||||
static integer i__, k;
|
||||
static doublereal q, q0;
|
||||
static integer ii;
|
||||
static doublereal qp;
|
||||
|
||||
/* Parameter adjustments */
|
||||
--l2;
|
||||
a_dim1 = *mp;
|
||||
a_offset = 1 + a_dim1;
|
||||
a -= a_offset;
|
||||
|
||||
/* Function Body */
|
||||
*ip = 0;
|
||||
if (*nl2 < 1) {
|
||||
return 0;
|
||||
}
|
||||
i__1 = *nl2;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
if (a[l2[i__] + 1 + (*kp + 1) * a_dim1] < -1e-6) {
|
||||
goto L2;
|
||||
}
|
||||
/* L11: */
|
||||
}
|
||||
return 0;
|
||||
L2:
|
||||
*q1 = -a[l2[i__] + 1 + a_dim1] / a[l2[i__] + 1 + (*kp + 1) * a_dim1];
|
||||
*ip = l2[i__];
|
||||
if (i__ + 1 > *nl2) {
|
||||
return 0;
|
||||
}
|
||||
i__1 = *nl2;
|
||||
for (++i__; i__ <= i__1; ++i__) {
|
||||
ii = l2[i__];
|
||||
if (a[ii + 1 + (*kp + 1) * a_dim1] < -1e-6) {
|
||||
q = -a[ii + 1 + a_dim1] / a[ii + 1 + (*kp + 1) * a_dim1];
|
||||
if (q < *q1) {
|
||||
*ip = ii;
|
||||
*q1 = q;
|
||||
} else if (q == *q1) {
|
||||
i__2 = *n;
|
||||
for (k = 1; k <= i__2; ++k) {
|
||||
qp = -a[*ip + 1 + (k + 1) * a_dim1] / a[*ip + 1 + (*kp +
|
||||
1) * a_dim1];
|
||||
q0 = -a[ii + 1 + (k + 1) * a_dim1] / a[ii + 1 + (*kp + 1)
|
||||
* a_dim1];
|
||||
if (q0 != qp) {
|
||||
goto L6;
|
||||
}
|
||||
/* L12: */
|
||||
}
|
||||
L6:
|
||||
if (q0 < qp) {
|
||||
*ip = ii;
|
||||
}
|
||||
}
|
||||
}
|
||||
/* L13: */
|
||||
}
|
||||
return 0;
|
||||
} /* simp2_ */
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,65 +0,0 @@
|
|||
/* simp3.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Subroutine */ int simp3_(doublereal *a, integer *mp, integer *np, integer *
|
||||
i1, integer *k1, integer *ip, integer *kp)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer a_dim1, a_offset, i__1, i__2;
|
||||
|
||||
/* Local variables */
|
||||
static integer ii, kk;
|
||||
static doublereal piv;
|
||||
|
||||
/* Parameter adjustments */
|
||||
a_dim1 = *mp;
|
||||
a_offset = 1 + a_dim1;
|
||||
a -= a_offset;
|
||||
|
||||
/* Function Body */
|
||||
piv = 1. / a[*ip + 1 + (*kp + 1) * a_dim1];
|
||||
if (*i1 >= 0) {
|
||||
i__1 = *i1 + 1;
|
||||
for (ii = 1; ii <= i__1; ++ii) {
|
||||
if (ii - 1 != *ip) {
|
||||
a[ii + (*kp + 1) * a_dim1] *= piv;
|
||||
i__2 = *k1 + 1;
|
||||
for (kk = 1; kk <= i__2; ++kk) {
|
||||
if (kk - 1 != *kp) {
|
||||
a[ii + kk * a_dim1] -= a[*ip + 1 + kk * a_dim1] * a[
|
||||
ii + (*kp + 1) * a_dim1];
|
||||
}
|
||||
/* L11: */
|
||||
}
|
||||
}
|
||||
/* L12: */
|
||||
}
|
||||
}
|
||||
i__1 = *k1 + 1;
|
||||
for (kk = 1; kk <= i__1; ++kk) {
|
||||
if (kk - 1 != *kp) {
|
||||
a[*ip + 1 + kk * a_dim1] = -a[*ip + 1 + kk * a_dim1] * piv;
|
||||
}
|
||||
/* L13: */
|
||||
}
|
||||
a[*ip + 1 + (*kp + 1) * a_dim1] = piv;
|
||||
return 0;
|
||||
} /* simp3_ */
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,199 +0,0 @@
|
|||
/* simplx.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
#include <stdio.h>
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static integer c__0 = 0;
|
||||
static integer c__1 = 1;
|
||||
|
||||
/* Subroutine */ int simplx_(doublereal *a, integer *m, integer *n, integer *
|
||||
mp, integer *np, integer *m1, integer *m2, integer *m3, integer *
|
||||
icase, integer *izrov, integer *iposv)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer a_dim1, a_offset, i__1, i__2;
|
||||
|
||||
/* Builtin functions */
|
||||
/* Subroutine */ int s_paus(char *, ftnlen);
|
||||
|
||||
/* Local variables */
|
||||
static integer i__, k, l1[1000], l2[1000], l3[1000];
|
||||
static doublereal q1;
|
||||
static integer m12, kh, ip, ir, kp, is, nl1, nl2;
|
||||
static doublereal bmax;
|
||||
extern /* Subroutine */ int simp1_(doublereal *, integer *, integer *,
|
||||
integer *, integer *, integer *, integer *, integer *, doublereal
|
||||
*), simp2_(doublereal *, integer *, integer *, integer *, integer
|
||||
*, integer *, integer *, integer *, integer *, doublereal *),
|
||||
simp3_(doublereal *, integer *, integer *, integer *, integer *,
|
||||
integer *, integer *);
|
||||
|
||||
/* Parameter adjustments */
|
||||
--iposv;
|
||||
--izrov;
|
||||
a_dim1 = *mp;
|
||||
a_offset = 1 + a_dim1;
|
||||
a -= a_offset;
|
||||
|
||||
/* Function Body */
|
||||
if (*m != *m1 + *m2 + *m3) {
|
||||
printf(" %d %d %d %d ", *m, *m1, *m2, *m3);
|
||||
s_paus("Bad input constraint counts.", (ftnlen)28);
|
||||
}
|
||||
nl1 = *n;
|
||||
i__1 = *n;
|
||||
for (k = 1; k <= i__1; ++k) {
|
||||
l1[k - 1] = k;
|
||||
izrov[k] = k;
|
||||
/* L11: */
|
||||
}
|
||||
nl2 = *m;
|
||||
i__1 = *m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
if (a[i__ + 1 + a_dim1] < 0.) {
|
||||
/* write(*,*) 'The A matrix input to SIMPLX is invalid.' */
|
||||
s_paus("Bad input tableau.", (ftnlen)18);
|
||||
}
|
||||
l2[i__ - 1] = i__;
|
||||
iposv[i__] = *n + i__;
|
||||
/* L12: */
|
||||
}
|
||||
i__1 = *m2;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
l3[i__ - 1] = 1;
|
||||
/* L13: */
|
||||
}
|
||||
ir = 0;
|
||||
if (*m2 + *m3 == 0) {
|
||||
goto L30;
|
||||
}
|
||||
ir = 1;
|
||||
i__1 = *n + 1;
|
||||
for (k = 1; k <= i__1; ++k) {
|
||||
q1 = (float)0.;
|
||||
i__2 = *m;
|
||||
for (i__ = *m1 + 1; i__ <= i__2; ++i__) {
|
||||
q1 += a[i__ + 1 + k * a_dim1];
|
||||
/* L14: */
|
||||
}
|
||||
a[*m + 2 + k * a_dim1] = -q1;
|
||||
/* L15: */
|
||||
}
|
||||
L10:
|
||||
i__1 = *m + 1;
|
||||
simp1_(&a[a_offset], mp, np, &i__1, l1, &nl1, &c__0, &kp, &bmax);
|
||||
if (bmax <= 1e-6 && a[*m + 2 + a_dim1] < -1e-6) {
|
||||
*icase = -1;
|
||||
return 0;
|
||||
} else if (bmax <= 1e-6 && a[*m + 2 + a_dim1] <= 1e-6) {
|
||||
m12 = *m1 + *m2 + 1;
|
||||
if (m12 <= *m) {
|
||||
i__1 = *m;
|
||||
for (ip = m12; ip <= i__1; ++ip) {
|
||||
if (iposv[ip] == ip + *n) {
|
||||
simp1_(&a[a_offset], mp, np, &ip, l1, &nl1, &c__1, &kp, &
|
||||
bmax);
|
||||
if (bmax > (float)0.) {
|
||||
goto L1;
|
||||
}
|
||||
}
|
||||
/* L16: */
|
||||
}
|
||||
}
|
||||
ir = 0;
|
||||
--m12;
|
||||
if (*m1 + 1 > m12) {
|
||||
goto L30;
|
||||
}
|
||||
i__1 = m12;
|
||||
for (i__ = *m1 + 1; i__ <= i__1; ++i__) {
|
||||
if (l3[i__ - *m1 - 1] == 1) {
|
||||
i__2 = *n + 1;
|
||||
for (k = 1; k <= i__2; ++k) {
|
||||
a[i__ + 1 + k * a_dim1] = -a[i__ + 1 + k * a_dim1];
|
||||
/* L17: */
|
||||
}
|
||||
}
|
||||
/* L18: */
|
||||
}
|
||||
goto L30;
|
||||
}
|
||||
simp2_(&a[a_offset], m, n, mp, np, l2, &nl2, &ip, &kp, &q1);
|
||||
if (ip == 0) {
|
||||
*icase = -1;
|
||||
return 0;
|
||||
}
|
||||
L1:
|
||||
i__1 = *m + 1;
|
||||
simp3_(&a[a_offset], mp, np, &i__1, n, &ip, &kp);
|
||||
if (iposv[ip] >= *n + *m1 + *m2 + 1) {
|
||||
i__1 = nl1;
|
||||
for (k = 1; k <= i__1; ++k) {
|
||||
if (l1[k - 1] == kp) {
|
||||
goto L2;
|
||||
}
|
||||
/* L19: */
|
||||
}
|
||||
L2:
|
||||
--nl1;
|
||||
i__1 = nl1;
|
||||
for (is = k; is <= i__1; ++is) {
|
||||
l1[is - 1] = l1[is];
|
||||
/* L21: */
|
||||
}
|
||||
} else {
|
||||
if (iposv[ip] < *n + *m1 + 1) {
|
||||
goto L20;
|
||||
}
|
||||
kh = iposv[ip] - *m1 - *n;
|
||||
if (l3[kh - 1] == 0) {
|
||||
goto L20;
|
||||
}
|
||||
l3[kh - 1] = 0;
|
||||
}
|
||||
a[*m + 2 + (kp + 1) * a_dim1] += (float)1.;
|
||||
i__1 = *m + 2;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
a[i__ + (kp + 1) * a_dim1] = -a[i__ + (kp + 1) * a_dim1];
|
||||
/* L22: */
|
||||
}
|
||||
L20:
|
||||
is = izrov[kp];
|
||||
izrov[kp] = iposv[ip];
|
||||
iposv[ip] = is;
|
||||
if (ir != 0) {
|
||||
goto L10;
|
||||
}
|
||||
L30:
|
||||
simp1_(&a[a_offset], mp, np, &c__0, l1, &nl1, &c__0, &kp, &bmax);
|
||||
if (bmax <= (float)0.) {
|
||||
*icase = 0;
|
||||
return 0;
|
||||
}
|
||||
simp2_(&a[a_offset], m, n, mp, np, l2, &nl2, &ip, &kp, &q1);
|
||||
if (ip == 0) {
|
||||
*icase = 1;
|
||||
return 0;
|
||||
}
|
||||
simp3_(&a[a_offset], mp, np, m, n, &ip, &kp);
|
||||
goto L20;
|
||||
} /* simplx_ */
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,65 +0,0 @@
|
|||
/* splie2.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static doublereal c_b4 = 1e30;
|
||||
|
||||
/* Subroutine */ int splie2_(doublereal *x1a, doublereal *x2a, doublereal *ya,
|
||||
integer *m, integer *n, doublereal *y2a)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer ya_dim1, ya_offset, y2a_dim1, y2a_offset, i__1, i__2;
|
||||
|
||||
/* Local variables */
|
||||
static integer j, k;
|
||||
static doublereal ytmp[100], y2tmp[100];
|
||||
extern /* Subroutine */ int spline_(doublereal *, doublereal *, integer *,
|
||||
doublereal *, doublereal *, doublereal *);
|
||||
|
||||
/* Parameter adjustments */
|
||||
--x1a;
|
||||
y2a_dim1 = *m;
|
||||
y2a_offset = 1 + y2a_dim1;
|
||||
y2a -= y2a_offset;
|
||||
ya_dim1 = *m;
|
||||
ya_offset = 1 + ya_dim1;
|
||||
ya -= ya_offset;
|
||||
--x2a;
|
||||
|
||||
/* Function Body */
|
||||
i__1 = *m;
|
||||
for (j = 1; j <= i__1; ++j) {
|
||||
i__2 = *n;
|
||||
for (k = 1; k <= i__2; ++k) {
|
||||
ytmp[k - 1] = ya[j + k * ya_dim1];
|
||||
/* L11: */
|
||||
}
|
||||
spline_(&x2a[1], ytmp, n, &c_b4, &c_b4, y2tmp);
|
||||
i__2 = *n;
|
||||
for (k = 1; k <= i__2; ++k) {
|
||||
y2a[j + k * y2a_dim1] = y2tmp[k - 1];
|
||||
/* L12: */
|
||||
}
|
||||
/* L13: */
|
||||
}
|
||||
return 0;
|
||||
} /* splie2_ */
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,66 +0,0 @@
|
|||
/* splin2.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static doublereal c_b4 = 1e30;
|
||||
|
||||
/* Subroutine */ int splin2_(doublereal *x1a, doublereal *x2a, doublereal *ya,
|
||||
doublereal *y2a, integer *m, integer *n, doublereal *x1, doublereal *
|
||||
x2, doublereal *y)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer ya_dim1, ya_offset, y2a_dim1, y2a_offset, i__1, i__2;
|
||||
|
||||
/* Local variables */
|
||||
static integer j, k;
|
||||
static doublereal ytmp[100], y2tmp[100], yytmp[100];
|
||||
extern /* Subroutine */ int spline_(doublereal *, doublereal *, integer *,
|
||||
doublereal *, doublereal *, doublereal *), splint_(doublereal *,
|
||||
doublereal *, doublereal *, integer *, doublereal *, doublereal *)
|
||||
;
|
||||
|
||||
/* Parameter adjustments */
|
||||
--x1a;
|
||||
y2a_dim1 = *m;
|
||||
y2a_offset = 1 + y2a_dim1;
|
||||
y2a -= y2a_offset;
|
||||
ya_dim1 = *m;
|
||||
ya_offset = 1 + ya_dim1;
|
||||
ya -= ya_offset;
|
||||
--x2a;
|
||||
|
||||
/* Function Body */
|
||||
i__1 = *m;
|
||||
for (j = 1; j <= i__1; ++j) {
|
||||
i__2 = *n;
|
||||
for (k = 1; k <= i__2; ++k) {
|
||||
ytmp[k - 1] = ya[j + k * ya_dim1];
|
||||
y2tmp[k - 1] = y2a[j + k * y2a_dim1];
|
||||
/* L11: */
|
||||
}
|
||||
splint_(&x2a[1], ytmp, y2tmp, n, x2, &yytmp[j - 1]);
|
||||
/* L12: */
|
||||
}
|
||||
spline_(&x1a[1], yytmp, m, &c_b4, &c_b4, y2tmp);
|
||||
splint_(&x1a[1], yytmp, y2tmp, m, x1, y);
|
||||
return 0;
|
||||
} /* splin2_ */
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,69 +0,0 @@
|
|||
/* spline.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Subroutine */ int spline_(doublereal *x, doublereal *y, integer *n,
|
||||
doublereal *yp1, doublereal *ypn, doublereal *y2)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer i__1;
|
||||
|
||||
/* Local variables */
|
||||
static integer i__, k;
|
||||
static doublereal p, u[100], qn, un, sig;
|
||||
|
||||
/* Parameter adjustments */
|
||||
--y2;
|
||||
--y;
|
||||
--x;
|
||||
|
||||
/* Function Body */
|
||||
if (*yp1 > 9.9e29) {
|
||||
y2[1] = (float)0.;
|
||||
u[0] = (float)0.;
|
||||
} else {
|
||||
y2[1] = -.5;
|
||||
u[0] = 3. / (x[2] - x[1]) * ((y[2] - y[1]) / (x[2] - x[1]) - *yp1);
|
||||
}
|
||||
i__1 = *n - 1;
|
||||
for (i__ = 2; i__ <= i__1; ++i__) {
|
||||
sig = (x[i__] - x[i__ - 1]) / (x[i__ + 1] - x[i__ - 1]);
|
||||
p = sig * y2[i__ - 1] + 2.;
|
||||
y2[i__] = (sig - (float)1.) / p;
|
||||
u[i__ - 1] = (((y[i__ + 1] - y[i__]) / (x[i__ + 1] - x[i__]) - (y[i__]
|
||||
- y[i__ - 1]) / (x[i__] - x[i__ - 1])) * 6. / (x[i__ + 1] -
|
||||
x[i__ - 1]) - sig * u[i__ - 2]) / p;
|
||||
/* L11: */
|
||||
}
|
||||
if (*ypn > 9.9e29) {
|
||||
qn = 0.;
|
||||
un = 0.;
|
||||
} else {
|
||||
qn = .5;
|
||||
un = 3. / (x[*n] - x[*n - 1]) * (*ypn - (y[*n] - y[*n - 1]) / (x[*n]
|
||||
- x[*n - 1]));
|
||||
}
|
||||
y2[*n] = (un - qn * u[*n - 2]) / (qn * y2[*n - 1] + 1.);
|
||||
for (k = *n - 1; k >= 1; --k) {
|
||||
y2[k] = y2[k] * y2[k + 1] + u[k - 1];
|
||||
/* L12: */
|
||||
}
|
||||
return 0;
|
||||
} /* spline_ */
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,68 +0,0 @@
|
|||
/* splint.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Subroutine */ int splint_(doublereal *xa, doublereal *ya, doublereal *y2a,
|
||||
integer *n, doublereal *x, doublereal *y)
|
||||
{
|
||||
/* System generated locals */
|
||||
doublereal d__1, d__2, d__3;
|
||||
|
||||
/* Builtin functions */
|
||||
/* Subroutine */ int s_paus(char *, ftnlen);
|
||||
|
||||
/* Local variables */
|
||||
static doublereal a, b, h__;
|
||||
static integer k, khi, klo;
|
||||
|
||||
/* Parameter adjustments */
|
||||
--y2a;
|
||||
--ya;
|
||||
--xa;
|
||||
|
||||
/* Function Body */
|
||||
klo = 1;
|
||||
khi = *n;
|
||||
L1:
|
||||
if (khi - klo > 1) {
|
||||
k = (khi + klo) / 2;
|
||||
if (xa[k] > *x) {
|
||||
khi = k;
|
||||
} else {
|
||||
klo = k;
|
||||
}
|
||||
goto L1;
|
||||
}
|
||||
h__ = xa[khi] - xa[klo];
|
||||
if (h__ == (float)0.) {
|
||||
s_paus("Bad XA input.", (ftnlen)13);
|
||||
}
|
||||
a = (xa[khi] - *x) / h__;
|
||||
b = (*x - xa[klo]) / h__;
|
||||
/* Computing 3rd power */
|
||||
d__1 = a;
|
||||
/* Computing 3rd power */
|
||||
d__2 = b;
|
||||
/* Computing 2nd power */
|
||||
d__3 = h__;
|
||||
*y = a * ya[klo] + b * ya[khi] + ((d__1 * (d__1 * d__1) - a) * y2a[klo] +
|
||||
(d__2 * (d__2 * d__2) - b) * y2a[khi]) * (d__3 * d__3) / 6.;
|
||||
return 0;
|
||||
} /* splint_ */
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,39 +0,0 @@
|
|||
# $License$
|
||||
#
|
||||
# $Id: Makefile.in,v 1.8 2008/12/30 21:49:42 hkmoffa Exp $
|
||||
#
|
||||
#/bin/sh
|
||||
.SUFFIXES :
|
||||
.SUFFIXES : .f .d .o
|
||||
|
||||
LIB = @buildlib@/librecipes.a
|
||||
|
||||
do_ranlib = @DO_RANLIB@
|
||||
|
||||
PURIFY=@PURIFY@
|
||||
|
||||
PIC_FLAG=@PIC@
|
||||
F_FLAGS = @FFLAGS@ $(PIC_FLAG)
|
||||
|
||||
OBJS = simp1.o simp2.o simp3.o simplx.o
|
||||
SRCS = $(OBJS:.o=.cpp)
|
||||
|
||||
|
||||
all: $(LIB)
|
||||
|
||||
$(LIB): $(OBJS)
|
||||
@ARCHIVE@ $(LIB) $(OBJS) > /dev/null
|
||||
ifeq ($(do_ranlib),1)
|
||||
@RANLIB@ $(LIB)
|
||||
endif
|
||||
|
||||
%.o : %.f
|
||||
$(PURIFY) @F77@ -c $< $(F_FLAGS)
|
||||
|
||||
clean:
|
||||
$(RM) $(OBJS) $(LIB)
|
||||
|
||||
depends:
|
||||
|
||||
|
||||
|
||||
|
|
@ -1,19 +0,0 @@
|
|||
SUBROUTINE SIMP1(A,MP,NP,MM,LL,NLL,IABF,KP,BMAX)
|
||||
implicit double precision (a-h,o-z)
|
||||
DIMENSION A(MP,NP),LL(NP)
|
||||
KP=LL(1)
|
||||
BMAX=A(MM+1,KP+1)
|
||||
IF(NLL.LT.2)RETURN
|
||||
DO 11 K=2,NLL
|
||||
IF(IABF.EQ.0)THEN
|
||||
TEST=A(MM+1,LL(K)+1)-BMAX
|
||||
ELSE
|
||||
TEST=ABS(A(MM+1,LL(K)+1))-ABS(BMAX)
|
||||
ENDIF
|
||||
IF(TEST.GT.0.)THEN
|
||||
BMAX=A(MM+1,LL(K)+1)
|
||||
KP=LL(K)
|
||||
ENDIF
|
||||
11 CONTINUE
|
||||
RETURN
|
||||
END
|
||||
|
|
@ -1,32 +0,0 @@
|
|||
SUBROUTINE SIMP2(A,M,N,MP,NP,L2,NL2,IP,KP,Q1)
|
||||
implicit double precision (a-h,o-z)
|
||||
PARAMETER (EPS=1.d-6)
|
||||
DIMENSION A(MP,NP),L2(MP)
|
||||
IP=0
|
||||
IF(NL2.LT.1)RETURN
|
||||
DO 11 I=1,NL2
|
||||
IF(A(L2(I)+1,KP+1).LT.-EPS)GO TO 2
|
||||
11 CONTINUE
|
||||
RETURN
|
||||
2 Q1=-A(L2(I)+1,1)/A(L2(I)+1,KP+1)
|
||||
IP=L2(I)
|
||||
IF(I+1.GT.NL2)RETURN
|
||||
DO 13 I=I+1,NL2
|
||||
II=L2(I)
|
||||
IF(A(II+1,KP+1).LT.-EPS)THEN
|
||||
Q=-A(II+1,1)/A(II+1,KP+1)
|
||||
IF(Q.LT.Q1)THEN
|
||||
IP=II
|
||||
Q1=Q
|
||||
ELSE IF (Q.EQ.Q1) THEN
|
||||
DO 12 K=1,N
|
||||
QP=-A(IP+1,K+1)/A(IP+1,KP+1)
|
||||
Q0=-A(II+1,K+1)/A(II+1,KP+1)
|
||||
IF(Q0.NE.QP)GO TO 6
|
||||
12 CONTINUE
|
||||
6 IF(Q0.LT.QP)IP=II
|
||||
ENDIF
|
||||
ENDIF
|
||||
13 CONTINUE
|
||||
RETURN
|
||||
END
|
||||
|
|
@ -1,22 +0,0 @@
|
|||
SUBROUTINE SIMP3(A,MP,NP,I1,K1,IP,KP)
|
||||
implicit double precision (a-h,o-z)
|
||||
DIMENSION A(MP,NP)
|
||||
PIV=1.d0/A(IP+1,KP+1)
|
||||
IF(I1.GE.0)THEN
|
||||
DO 12 II=1,I1+1
|
||||
IF(II-1.NE.IP)THEN
|
||||
A(II,KP+1)=A(II,KP+1)*PIV
|
||||
DO 11 KK=1,K1+1
|
||||
IF(KK-1.NE.KP)THEN
|
||||
A(II,KK)=A(II,KK)-A(IP+1,KK)*A(II,KP+1)
|
||||
ENDIF
|
||||
11 CONTINUE
|
||||
ENDIF
|
||||
12 CONTINUE
|
||||
ENDIF
|
||||
DO 13 KK=1,K1+1
|
||||
IF(KK-1.NE.KP)A(IP+1,KK)=-A(IP+1,KK)*PIV
|
||||
13 CONTINUE
|
||||
A(IP+1,KP+1)=PIV
|
||||
RETURN
|
||||
END
|
||||
|
|
@ -1,99 +0,0 @@
|
|||
SUBROUTINE SIMPLX(A,M,N,MP,NP,M1,M2,M3,ICASE,IZROV,IPOSV)
|
||||
implicit double precision (a-h,o-z)
|
||||
PARAMETER(MMAX=1000,EPS=1.d-6)
|
||||
DIMENSION A(MP,NP),IZROV(N),IPOSV(M),L1(MMAX),L2(MMAX),L3(MMAX)
|
||||
IF(M.NE.M1+M2+M3)PAUSE 'Bad input constraint counts.'
|
||||
NL1=N
|
||||
DO 11 K=1,N
|
||||
L1(K)=K
|
||||
IZROV(K)=K
|
||||
11 CONTINUE
|
||||
NL2=M
|
||||
DO 12 I=1,M
|
||||
IF(A(I+1,1).LT.0.d0) then
|
||||
c write(*,*) 'The A matrix input to SIMPLX is invalid.'
|
||||
PAUSE 'Bad input tableau.'
|
||||
end if
|
||||
L2(I)=I
|
||||
IPOSV(I)=N+I
|
||||
12 CONTINUE
|
||||
DO 13 I=1,M2
|
||||
L3(I)=1
|
||||
13 CONTINUE
|
||||
IR=0
|
||||
IF(M2+M3.EQ.0)GO TO 30
|
||||
IR=1
|
||||
DO 15 K=1,N+1
|
||||
Q1=0.
|
||||
DO 14 I=M1+1,M
|
||||
Q1=Q1+A(I+1,K)
|
||||
14 CONTINUE
|
||||
A(M+2,K)=-Q1
|
||||
15 CONTINUE
|
||||
10 CALL SIMP1(A,MP,NP,M+1,L1,NL1,0,KP,BMAX)
|
||||
IF(BMAX.LE.EPS.AND.A(M+2,1).LT.-EPS)THEN
|
||||
ICASE=-1
|
||||
RETURN
|
||||
ELSE IF(BMAX.LE.EPS.AND.A(M+2,1).LE.EPS)THEN
|
||||
M12=M1+M2+1
|
||||
IF(M12.LE.M)THEN
|
||||
DO 16 IP=M12,M
|
||||
IF(IPOSV(IP).EQ.IP+N)THEN
|
||||
CALL SIMP1(A,MP,NP,IP,L1,NL1,1,KP,BMAX)
|
||||
IF(BMAX.GT.0.)GO TO 1
|
||||
ENDIF
|
||||
16 CONTINUE
|
||||
ENDIF
|
||||
IR=0
|
||||
M12=M12-1
|
||||
IF(M1+1.GT.M12)GO TO 30
|
||||
DO 18 I=M1+1,M12
|
||||
IF(L3(I-M1).EQ.1)THEN
|
||||
DO 17 K=1,N+1
|
||||
A(I+1,K)=-A(I+1,K)
|
||||
17 CONTINUE
|
||||
ENDIF
|
||||
18 CONTINUE
|
||||
GO TO 30
|
||||
ENDIF
|
||||
CALL SIMP2(A,M,N,MP,NP,L2,NL2,IP,KP,Q1)
|
||||
IF(IP.EQ.0)THEN
|
||||
ICASE=-1
|
||||
RETURN
|
||||
ENDIF
|
||||
1 CALL SIMP3(A,MP,NP,M+1,N,IP,KP)
|
||||
IF(IPOSV(IP).GE.N+M1+M2+1)THEN
|
||||
DO 19 K=1,NL1
|
||||
IF(L1(K).EQ.KP)GO TO 2
|
||||
19 CONTINUE
|
||||
2 NL1=NL1-1
|
||||
DO 21 IS=K,NL1
|
||||
L1(IS)=L1(IS+1)
|
||||
21 CONTINUE
|
||||
ELSE
|
||||
IF(IPOSV(IP).LT.N+M1+1)GO TO 20
|
||||
KH=IPOSV(IP)-M1-N
|
||||
IF(L3(KH).EQ.0)GO TO 20
|
||||
L3(KH)=0
|
||||
ENDIF
|
||||
A(M+2,KP+1)=A(M+2,KP+1)+1.
|
||||
DO 22 I=1,M+2
|
||||
A(I,KP+1)=-A(I,KP+1)
|
||||
22 CONTINUE
|
||||
20 IS=IZROV(KP)
|
||||
IZROV(KP)=IPOSV(IP)
|
||||
IPOSV(IP)=IS
|
||||
IF(IR.NE.0)GO TO 10
|
||||
30 CALL SIMP1(A,MP,NP,0,L1,NL1,0,KP,BMAX)
|
||||
IF(BMAX.LE.0.)THEN
|
||||
ICASE=0
|
||||
RETURN
|
||||
ENDIF
|
||||
CALL SIMP2(A,M,N,MP,NP,L2,NL2,IP,KP,Q1)
|
||||
IF(IP.EQ.0)THEN
|
||||
ICASE=1
|
||||
RETURN
|
||||
ENDIF
|
||||
CALL SIMP3(A,MP,NP,M,N,IP,KP)
|
||||
GO TO 20
|
||||
END
|
||||
|
|
@ -1,15 +0,0 @@
|
|||
SUBROUTINE SPLIE2(X1A,X2A,YA,M,N,Y2A)
|
||||
implicit double precision (a-h,o-z)
|
||||
PARAMETER (NN=100)
|
||||
DIMENSION X1A(M),X2A(N),YA(M,N),Y2A(M,N),YTMP(NN),Y2TMP(NN)
|
||||
DO 13 J=1,M
|
||||
DO 11 K=1,N
|
||||
YTMP(K)=YA(J,K)
|
||||
11 CONTINUE
|
||||
CALL SPLINE(X2A,YTMP,N,1.D30,1.D30,Y2TMP)
|
||||
DO 12 K=1,N
|
||||
Y2A(J,K)=Y2TMP(K)
|
||||
12 CONTINUE
|
||||
13 CONTINUE
|
||||
RETURN
|
||||
END
|
||||
|
|
@ -1,16 +0,0 @@
|
|||
SUBROUTINE SPLIN2(X1A,X2A,YA,Y2A,M,N,X1,X2,Y)
|
||||
implicit double precision (a-h,o-z)
|
||||
PARAMETER (NN=100)
|
||||
DIMENSION X1A(M),X2A(N),YA(M,N),Y2A(M,N),YTMP(NN),
|
||||
$ Y2TMP(NN),YYTMP(NN)
|
||||
DO 12 J=1,M
|
||||
DO 11 K=1,N
|
||||
YTMP(K)=YA(J,K)
|
||||
Y2TMP(K)=Y2A(J,K)
|
||||
11 CONTINUE
|
||||
CALL SPLINT(X2A,YTMP,Y2TMP,N,X2,YYTMP(J))
|
||||
12 CONTINUE
|
||||
CALL SPLINE(X1A,YYTMP,M,1.d30,1.d30,Y2TMP)
|
||||
CALL SPLINT(X1A,YYTMP,Y2TMP,M,X1,Y)
|
||||
RETURN
|
||||
END
|
||||
|
|
@ -1,31 +0,0 @@
|
|||
SUBROUTINE SPLINE(X,Y,N,YP1,YPN,Y2)
|
||||
implicit double precision (a-h,o-z)
|
||||
PARAMETER (NMAX=100)
|
||||
DIMENSION X(N),Y(N),Y2(N),U(NMAX)
|
||||
IF (YP1.GT..99D30) THEN
|
||||
Y2(1)=0.
|
||||
U(1)=0.
|
||||
ELSE
|
||||
Y2(1)=-0.5d0
|
||||
U(1)=(3.d0/(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1)
|
||||
ENDIF
|
||||
DO 11 I=2,N-1
|
||||
SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1))
|
||||
P=SIG*Y2(I-1)+2.d0
|
||||
Y2(I)=(SIG-1.)/P
|
||||
U(I)=(6.d0*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))
|
||||
* /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P
|
||||
11 CONTINUE
|
||||
IF (YPN.GT..99d30) THEN
|
||||
QN=0.d0
|
||||
UN=0.d0
|
||||
ELSE
|
||||
QN=0.5d0
|
||||
UN=(3.d0/(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1)))
|
||||
ENDIF
|
||||
Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.d0)
|
||||
DO 12 K=N-1,1,-1
|
||||
Y2(K)=Y2(K)*Y2(K+1)+U(K)
|
||||
12 CONTINUE
|
||||
RETURN
|
||||
END
|
||||
|
|
@ -1,22 +0,0 @@
|
|||
SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y)
|
||||
implicit double precision (a-h,o-z)
|
||||
DIMENSION XA(N),YA(N),Y2A(N)
|
||||
KLO=1
|
||||
KHI=N
|
||||
1 IF (KHI-KLO.GT.1) THEN
|
||||
K=(KHI+KLO)/2
|
||||
IF(XA(K).GT.X)THEN
|
||||
KHI=K
|
||||
ELSE
|
||||
KLO=K
|
||||
ENDIF
|
||||
GOTO 1
|
||||
ENDIF
|
||||
H=XA(KHI)-XA(KLO)
|
||||
IF (H.EQ.0.) PAUSE 'Bad XA input.'
|
||||
A=(XA(KHI)-X)/H
|
||||
B=(X-XA(KLO))/H
|
||||
Y=A*YA(KLO)+B*YA(KHI)+
|
||||
* ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6.d0
|
||||
RETURN
|
||||
END
|
||||
|
|
@ -1,34 +0,0 @@
|
|||
// Ideal Gas
|
||||
|
||||
#include "ideal.h"
|
||||
#include <math.h>
|
||||
|
||||
double ideal::up(){
|
||||
|
||||
}
|
||||
|
||||
double ideal::sp() {
|
||||
|
||||
|
||||
double ideal::Pp(){
|
||||
|
||||
|
||||
//equation s3
|
||||
double ideal::Psat(){
|
||||
|
||||
}
|
||||
|
||||
|
||||
//equation D3
|
||||
double ideal::ldens(){
|
||||
|
||||
}
|
||||
|
||||
double methane::Tcrit() {return 0.0;}
|
||||
double methane::Pcrit() {return 0.0;}
|
||||
double methane::Vcrit() {return 0.0;}
|
||||
double methane::Tmin() {return Tmn;}
|
||||
double methane::Tmax() {return Tmx;}
|
||||
char * methane::name() {return "ideal gas";}
|
||||
char * methane::formula() {return "--";}
|
||||
double methane::MolWt() {return M;}
|
||||
|
|
@ -1,80 +0,0 @@
|
|||
// Lee-Kesler equation of state
|
||||
|
||||
#include "lkw.h"
|
||||
#include <math.h>
|
||||
|
||||
const double omega_ref = 0.3978;
|
||||
|
||||
//--------------------------- member functions ------------------
|
||||
|
||||
double lkw::z(double Temp, double Pres) {
|
||||
T = Temp;
|
||||
f0.Set(TP,Temp,Pres);
|
||||
f1.Set(TP,Temp,Pres);
|
||||
if (Temp < Tcrit() && f0.x() != f1.x()) { // need to find Psat
|
||||
if (Pres >= Ps()) { //liquid
|
||||
f0.Set_meta(Liquid,Pres);
|
||||
f1.Set_meta(Liquid,Pres);
|
||||
}
|
||||
else {
|
||||
f0.Set_meta(Vapor,Pres);
|
||||
f1.Set_meta(Vapor,Pres);
|
||||
}
|
||||
}
|
||||
double z0 = f0.z();
|
||||
double zz = z0 + (f1.z() - z0)*Acent/omega_ref;
|
||||
return zz;
|
||||
}
|
||||
|
||||
|
||||
double lkw::Ps() {
|
||||
|
||||
// start with linear estimate of Psat
|
||||
double p0 = f0.Psat();
|
||||
double pse = p0 + (f1.Psat() - p0)/omega_ref;
|
||||
double lps = log(pse);
|
||||
double zf, zv;
|
||||
// loop until g_f = g_v
|
||||
for (int i = 0; i<20; i++) {
|
||||
|
||||
f0.Set_meta(Liquid,pse); // set both to liquid
|
||||
f1.Set_meta(Liquid,pse);
|
||||
double gf = f0.gp() + (f1.gp() - f0.gp())/omega_ref;
|
||||
zf = f0.z() + (f1.z() - f0.z())/omega_ref;
|
||||
|
||||
f0.Set_meta(Vapor,pse); // set both to vapor
|
||||
f1.Set_meta(Vapor,pse);
|
||||
double gv = f0.gp() + (f1.gp() - f0.gp())/omega_ref;
|
||||
zv = f0.z() + (f1.z() - f0.z())/omega_ref;
|
||||
|
||||
double gfv = gv - gf;
|
||||
double dlp = gfv*Mw/(8314.3*T*(zv - zf));
|
||||
lps -= dlp;
|
||||
pse = exp(lps);
|
||||
if (fabs(gfv) < 0.001) break;
|
||||
}
|
||||
if (i >= 20) {
|
||||
Pst = Undefined;
|
||||
Rhv = Undefined;
|
||||
Rhf = Undefined;
|
||||
Tslast = Undefined;
|
||||
set_Err(NoConverge);
|
||||
}
|
||||
else {
|
||||
Pst = pse;
|
||||
Tslast = T;
|
||||
Rhv = pse*Mw/(zv*8314.3*T);
|
||||
Rhf = pse*Mw/(zf*8314.3*T);
|
||||
}
|
||||
return Pst;
|
||||
}
|
||||
/*
|
||||
double lk::Tcrit() {return Tcr;}
|
||||
double lk::Pcrit() {return Pcr;}
|
||||
double lk::Vcrit() {return 0.2901*8314.3*Tcr/(Pcr*Mw);}
|
||||
double lk::Tmin() {return -100.0;}
|
||||
double lk::Tmax() {return 10000.0;}
|
||||
char * lk::name() {return "Lee-Kesler";}
|
||||
char * lk::formula() {return "---";}
|
||||
double lk::MolWt() {return Mw;}
|
||||
*/
|
||||
|
|
@ -1,49 +0,0 @@
|
|||
#ifndef LKW_H
|
||||
#define LKW_H
|
||||
|
||||
#include "lk.h"
|
||||
|
||||
class lkw {
|
||||
public:
|
||||
lkw(double tc = 1.0, double pc = 1.0, double wt = 8314.3, double omega = 0.0)
|
||||
{
|
||||
Tcr = tc;
|
||||
Pcr = pc;
|
||||
Mw = wt;
|
||||
Acent = omega;
|
||||
T = Undefined;
|
||||
Pressure = Undefined;
|
||||
};
|
||||
~lkw() {};
|
||||
/*
|
||||
double MolWt();
|
||||
double Tcrit();
|
||||
double Pcrit();
|
||||
double Vcrit();
|
||||
double Tmin();
|
||||
double Tmax();
|
||||
char * name();
|
||||
char * formula();
|
||||
*/
|
||||
double Ps();
|
||||
// double Pp();
|
||||
// double up();
|
||||
// double sp();
|
||||
// double Psat();
|
||||
double z(double tstar, double pstar);
|
||||
// double hdep();
|
||||
// double sdep();
|
||||
// double ldens();
|
||||
|
||||
protected:
|
||||
double Tcr, Pcr, Mw;
|
||||
double Acent;
|
||||
double Pressure;
|
||||
double T, Rho, Rhov, Rhof;
|
||||
lk f0; //(1.0, 1.0, 8314.3, 0);
|
||||
lk f1; //(1.0, 1.0, 8314.3, 1);
|
||||
|
||||
private:
|
||||
|
||||
};
|
||||
#endif // ! LKW_H
|
||||
Loading…
Add table
Reference in a new issue