home *** CD-ROM | disk | FTP | other *** search
/ ftp.ee.lbl.gov / 2014.05.ftp.ee.lbl.gov.tar / ftp.ee.lbl.gov / bmd-1.0beta.tar.Z / bmd-1.0beta.tar / bmd-1.0beta / app / omtd / pt.c < prev    next >
C/C++ Source or Header  |  1991-01-13  |  3KB  |  135 lines
  1. /*
  2.  * Copyright (c) 1990 Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms are permitted
  6.  * provided that the above copyright notice and this paragraph are
  7.  * duplicated in all such forms and that any documentation,
  8.  * advertising materials, and other materials related to such
  9.  * distribution and use acknowledge that the software was developed
  10.  * by the University of California, Lawrence Berkeley Laboratory,
  11.  * Berkeley, CA.  The name of the University may not be used to
  12.  * endorse or promote products derived from this software without
  13.  * specific prior written permission.
  14.  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
  15.  * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
  16.  * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
  17.  */
  18. /*
  19.  * Piece-wise continuous functions.
  20.  */
  21. #include <sys/types.h>
  22.  
  23. struct pt_list {
  24.     u_long x;
  25.     float y;
  26.     struct pt_list *next;
  27. };
  28.  
  29. pt_insert(p, x, y)
  30.     struct pt_list **p;
  31.     u_long x;
  32.     float y;
  33. {
  34.     struct pt_list *n = (struct pt_list *)malloc(sizeof(*n));
  35.  
  36.     n->x = x;
  37.     n->y = y;
  38.  
  39.     while (*p && x > (*p)->x)
  40.         p = &(*p)->next;
  41.     if (*p && (*p)->x == x) {
  42.         free((char *)n);
  43.         return;
  44.     }
  45.     n->next = *p;
  46.     *p = n;
  47. }
  48.  
  49. u_long
  50. pt_integrate(p, x0)
  51.     struct pt_list *p;
  52.     u_long x0;
  53. {
  54.     struct pt_list *n;
  55.     float sum = 0.0, h;
  56.  
  57.     while ((n = p->next) && x0 > n->x) {
  58.         sum += (p->y + n->y) * (n->x - p->x) / 2.0;
  59.         p = n;
  60.     }
  61.     if (n == 0 || p->y == n->y)
  62.         return sum + p->y * (x0 - p->x);
  63.     else {
  64.         h = p->y + (x0 - p->x) * (n->y - p->y) / (n->x - p->x);
  65.         return sum +  0.5 * (p->y + h) * (x0 - p->x);
  66.     }
  67. }
  68.  
  69. double sqrt();
  70. #define SQRT(x) (sqrt((double)(x)))
  71.  
  72. u_long
  73. pt_back_integrate(p, A)
  74.     struct pt_list *p;
  75.     u_long A;
  76. {
  77.     struct pt_list *n;
  78.     float sum = A, a, b;
  79.  
  80.     while (n = p->next) {
  81.         a = (p->y + n->y) * (n->x - p->x) / 2.0;
  82.         if (sum < a)
  83.             break;
  84.         sum -= a;
  85.         p = n;
  86.     }
  87.     if (n == 0 || p->y == n->y)
  88.         return p->x + (sum / p->y);
  89.     /* Solve quadratic. */
  90.     a = 0.5 * (n->y - p->y) / (n->x - p->x);
  91.     b = p->y;
  92.     return p->x + (SQRT(b * b + 4.0 * a * sum) - b) * 0.5 / a;
  93. }
  94.  
  95. pt_free(p)
  96.     struct pt_list *p;
  97. {
  98.     while (p) {
  99.         struct pt_list *next = p->next;
  100.         free((char *)p);
  101.         p = next;
  102.     }
  103. }
  104.  
  105. write_ptlist(fd, list)
  106.     int fd;
  107.     struct pt_list *list;
  108. {
  109.     struct pt_list *p;
  110.     int cnt = 0;
  111.  
  112.     for (p = list; p != 0; p = p->next)
  113.         ++cnt;
  114.     (void)write(fd, &cnt, sizeof cnt);
  115.     for (p = list; p != 0; p = p->next)
  116.         (void)write(fd, (char *)p, sizeof *p - sizeof p->next);
  117. }
  118.  
  119. read_ptlist(fd, listp)
  120.     int fd;
  121.     struct pt_list **listp;
  122. {
  123.     struct pt_list *p;
  124.     int cnt = 0;
  125.  
  126.     (void)read(fd, &cnt, sizeof cnt);
  127.     while (--cnt >= 0) {
  128.         p = (struct pt_list *)malloc(sizeof(*p));
  129.         (void)read(fd, (char *)p, sizeof *p - sizeof p->next);
  130.         *listp = p;
  131.         listp = &p->next;
  132.     }
  133.     *listp = 0;
  134. }
  135.