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