/*****************************************************************************

* .FILE  : numeric.c

* .TITLE : Расчет двойного интеграла при помощи метода Симпсона (парабол)

* .DESCR : Курсовой проект по численным методам (1994 год)

*        :

*        :  Выполнил: Петренко В.С. (гр. ПС-301)

*        :  Проверил: Панюков А.В.

*        :

* .NOTE  : NOT FOR RENTAL OR SALE.

*        : FEDERAL LAW PROVIDES SEVERE CIVIL & CRIMINAL PENALTIES FOR

*        : UNAUTHORIZED DUPLICATION OR DISTRIBUTION.

*        :

*        : (C) '94 by P$P

*****************************************************************************/

#include <math.h>

#include <stdio.h>

/*****************************************************************************

* .NAME   : m_Simpson

* .TITLE  : Расчет интеграла методом Симпсона (парабол)

* .DESCR  :

*         :

* .PARAMS : double m_Simpson (double (*func) (double, double),

*         :                   double t_fix, double t_limit, int N);

*         :        double (*func) (double, double) - подынтегральная ф-я

*         :        double t_fix   - фиксированный первый аргумент

*         :        double t_limit - верхний предел интегрирования,

*         :                         нижний равен -t_limit

*         :        int N          - число точек разбиения

* .RETURN : Значение вычисленного интеграла

*****************************************************************************/

double m_Simpson (double (*func) (double, double),

                  double t_fix, double t_limit, int N)

 {

   double sum1 = 0; /* -¬                        */

   double sum2 = 0; /*  ¦ временные переменные   */

   double sum3 = 0; /* --                        */

   double sum;      /* конечный результат        */

   double h    = (2 * t_limit) / N; /* шаг сетки */

   int i;           /* временная                 */
    sum1 = (*func) (t_fix, -t_limit) + (*func) (t_fix, +t_limit);
    for (i = 1; i <= N - 1; i++) sum2 += (*func) (t_fix, -t_limit + (i * h));

    sum2 *= 2;
    for (i = 1; i <= N; i++) sum3 += (*func) (t_fix, -t_limit + ((i-0.5) * h));

    sum3 *= 4;
    sum = sum1 + sum2 + sum3;

    sum = (h / 6) * sum;

  return sum;

 }

                      /*                       */

                      /* Глобальные переменные */

                      /* ~~~~~~~~~~~~~~~~~~~~~ */

#define PI 3.1415926536    /* число П                                   */

double k;                  /* параметр функции - задается пользователем */

int N_MAX;                 /* число узлов сетки разбиения     */

double (*currFunc) (double, double); /* выбранная пользователем функция */
double f1 (double x, double y);

double f2 (double x, double y);

double f3 (double x, double y);
/*****************************************************************************

* .NAME   : double F (double dummy, double t)

* .TITLE  : Вычисляет внутренний интеграл (G (t)).

*         : См. текст курсового проекта.

* .DESCR  :                           П/2

*         :                            -

*         : первая вычисляемая функция ¦ G (t + П/2) * sin (t + П/2) dt

*         :             t              -

*         :             -            -П/2

*         : где G (t) = ¦ currFunc (t, tau) dtau

*         :             -

*         :            -t

* .PARAMS : double F (double dummy, double t);

*         :        double dummy - фиктивный первый аргумент, при вызове этой

*         :                       функции он не используется, т.к. она

*         :                       функция одного аргумента

*         :        double       - действительный второй аргумент

* .RETURN : Значение функции: G (t) * sin (k * t);

*****************************************************************************/

double F (double dummy, double t)

 {

   double G;

   t = t + PI / 2; /* сдвижка начала координат, чтобы пределы    */

                   /* были симметричны (в нашем случае - на П/2) */

    G = m_Simpson (currFunc, t, t, N_MAX);

  return G * sin (k * t);

 }

/*****************************************************************************

* .NAME   : main

* .TITLE  : Основная диалоговая функция.

* .DESCR  : Запрашивается интересующая пользователя функция,

*         : параметр k и число узлов сетки N_MAX.

*         : Выводит на экран вычисленное значение интеграла и

*         : два справочных значения - П и П/2.

*         :

* .PARAMS : void main (void);

* .RETURN :

*****************************************************************************/

void main (void)

 {

   double integral;  /* значение вычисленного интеграла */

   int    selection; /* номер выбранной функции         */

   /* массив доступных функций */

   double (*functions []) (double, double) = { f1, f2, f3 };
    printf ("\n   Вычисление интеграла методом Симпсона (парабол)   ");

    printf ("\n   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   ");

    printf ("\n              --                                     ");

    printf ("\n          I = ¦¦ sin k(x + y) f (x, y) dx dy         ");

    printf ("\n              --                                     ");

    printf ("\n              D                                      ");

    printf ("\n где D = { (x, y): x, y >= 0; x + y <= П }, f Е C (D)");

    printf ("\n");

    printf ("\nДля какой функции рассчитывать:             ");

    printf ("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~             ");

    printf ("\n  1) f (x, y) = 0.5 * cos (y)               ");

    printf ("\n     -                        -- 0; k != 1  ");

    printf ("\n     ¦ sin x * sin (kx) dx => ¦             ");

    printf ("\n     -                        L- П/2; k = 1 ");

    printf ("\n                                            ");

    printf ("\n  2) f (x, y) = 0.5 - sin (y)               ");

    printf ("\n     -                                      ");

    printf ("\n     ¦ x * sin (kx) dx =====> П; k = 1      ");

    printf ("\n     -                                      ");

    printf ("\n                                            ");

    printf ("\n  3) f (x, y) = sqrt (x * x + y * y)");

    printf ("\n");

    do

     {

       printf ("Ваш выбор: ");

       scanf ("%d", &selection);

     } while (!(1 <= selection && selection <= 3));
    printf ("Параметр k: ");

    scanf ("%lg", &k);
    do

     {

       printf ("Число узлов сетки N: ");

       scanf ("%d", &N_MAX);

     } while (!(N_MAX > 0));

    printf ("\n");

    printf ("\n Расчет интеграла ...");
    currFunc = functions [selection - 1];          /* текущая функция    */

    integral = m_Simpson (F, 0, PI / 2, N_MAX);    /* вычисляем интеграл */

    printf ("\n Значение интеграла равно: %.12lg", integral); /* вывод   */

    printf ("\n Величины: П = %.12lg; П/2 = %.12lg", PI, PI / 2);

 }
/*****************************************************************************

* .FILE  : func.c

* .TITLE : Содержит функции пользователя, которые можно

*        : изменять без перекомпиляции основной программы

* .DESCR : После изменения этого модуля его необходимо перекомпилировать

*        : и слинковать с numeric.obj

*        :

*        :

* .NOTE  : NOT FOR RENTAL OR SALE.

*        : FEDERAL LAW PROVIDES SEVERE CIVIL & CRIMINAL PENALTIES FOR

*        : UNAUTHORIZED DUPLICATION OR DISTRIBUTION.

*        :

*        : (C) '94 by P$P

*****************************************************************************/

#include <math.h>
/* выбираемая пользователем функция No.1 */

double f1 (double x, double y)

 { return 0.5 * cos (y); }
/* выбираемая пользователем функция No.2 */

double f2 (double x, double y)

 { return 0.5 - sin (y); }
/* выбираемая пользователем функция No.3 */

double f3 (double x, double y)

 { return sqrt (x * x + y * y); }