single-laplace.c
1 //! @file single-laplace.c
2 //! @author J. Marcel van der Veer
3 //!
4 //! @section Copyright
5 //!
6 //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
7 //! Copyright 2001-2023 J. Marcel van der Veer [algol68g@xs4all.nl].
8 //!
9 //! @section License
10 //!
11 //! This program is free software; you can redistribute it and/or modify it
12 //! under the terms of the GNU General Public License as published by the
13 //! Free Software Foundation; either version 3 of the License, or
14 //! (at your option) any later version.
15 //!
16 //! This program is distributed in the hope that it will be useful, but
17 //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
19 //! more details. You should have received a copy of the GNU General Public
20 //! License along with this program. If not, see [http://www.gnu.org/licenses/].
21
22 //! @section Synopsis
23 //!
24 //! REAL laplace routines.
25
26 #include "a68g.h"
27 #include "a68g-genie.h"
28 #include "a68g-prelude.h"
29
30 #if defined (HAVE_GSL)
31
32 //! @brief Map GSL error handler onto a68g error handler.
33
34 void laplace_error_handler (const char *reason, const char *file, int line, int gsl_errno)
35 {
36 if (line != 0) {
37 ASSERT (snprintf (A68 (edit_line), SNPRINTF_SIZE, "%s in line %d of file %s", reason, line, file) >= 0);
38 } else {
39 ASSERT (snprintf (A68 (edit_line), SNPRINTF_SIZE, "%s", reason) >= 0);
40 }
41 diagnostic (A68_RUNTIME_ERROR, A68 (f_entry), ERROR_LAPLACE, A68 (edit_line), gsl_strerror (gsl_errno));
42 exit_genie (A68 (f_entry), A68_RUNTIME_ERROR);
43 }
44
45 //! @brief Detect math errors.
46
47 void laplace_test_error (int rc)
48 {
49 if (rc != 0) {
50 laplace_error_handler ("math error", "", 0, rc);
51 }
52 }
53
54 //! @brief PROC (PROC (REAL) REAL, REAL, REF REAL) REAL laplace
55
56 #define LAPLACE_DIVISIONS 1024
57
58 typedef struct A68_LAPLACE A68_LAPLACE;
59
60 struct A68_LAPLACE
61 {
62 NODE_T *p;
63 A68_PROCEDURE f;
64 REAL_T s;
65 };
66
67 //! @brief Evaluate function for Laplace transform.
68
69 REAL_T laplace_f (REAL_T t, void *z)
70 {
71 A68_LAPLACE *l = (A68_LAPLACE *) z;
72 ADDR_T pop_sp = A68_SP, pop_fp = A68_FP;
73 MOID_T *u = M_PROC_REAL_REAL;
74 A68_REAL *ft = (A68_REAL *) STACK_TOP;
75 PUSH_VALUE (P (l), t, A68_REAL);
76 genie_call_procedure (P (l), MOID (&(F (l))), u, u, &(F (l)), pop_sp, pop_fp);
77 A68_SP = pop_sp;
78 return VALUE (ft) * a68_exp (-(S (l)) * t);
79 }
80
81 //! @brief Calculate Laplace transform.
82
83 void genie_laplace (NODE_T * p)
84 {
85 A68_REF ref_error;
86 A68_REAL s, *error;
87 A68_PROCEDURE f;
88 A68_LAPLACE l;
89 gsl_function g;
90 gsl_integration_workspace *w;
91 REAL_T result, estimated_error;
92 int rc;
93 gsl_error_handler_t *save_handler = gsl_set_error_handler (laplace_error_handler);
94 POP_REF (p, &ref_error);
95 CHECK_REF (p, ref_error, M_REF_REAL);
96 error = (A68_REAL *) ADDRESS (&ref_error);
97 POP_OBJECT (p, &s, A68_REAL);
98 POP_PROCEDURE (p, &f);
99 P (&l) = p;
100 F (&l) = f;
101 S (&l) = VALUE (&s);
102 FUNCTION (&g) = &laplace_f;
103 GSL_PARAMS (&g) = &l;
104 w = gsl_integration_workspace_alloc (LAPLACE_DIVISIONS);
105 if (VALUE (error) >= 0.0) {
106 rc = gsl_integration_qagiu (&g, 0.0, VALUE (error), 0.0, LAPLACE_DIVISIONS, w, &result, &estimated_error);
107 } else {
108 rc = gsl_integration_qagiu (&g, 0.0, 0.0, -VALUE (error), LAPLACE_DIVISIONS, w, &result, &estimated_error);
109 }
110 laplace_test_error (rc);
111 VALUE (error) = estimated_error;
112 PUSH_VALUE (p, result, A68_REAL);
113 gsl_integration_workspace_free (w);
114 (void) gsl_set_error_handler (save_handler);
115 }
116
117 #endif
