File: derivative.c

package info (click to toggle)
xgraph 12.1-13
  • links: PTS, VCS
  • area: main
  • in suites: squeeze
  • size: 1,472 kB
  • ctags: 1,937
  • sloc: ansic: 7,924; sh: 3,464; makefile: 55
file content (130 lines) | stat: -rw-r--r-- 4,362 bytes parent folder | download | duplicates (5)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
/* This entire routine written by PW. */

#include <stdio.h>
#include <string.h>
#include "xgraph.h"

void
Der1()
{
  PointList *theList, *D1List, *D2List;
  double m,b;
  int i;

  
  for (i=0;i<MAXSETS;i++) {
    theList = PlotData[i].list;
    if (theList) {
      DataD1[i].list = (PointList *)malloc(sizeof(PointList));
      DataD2[i].list = (PointList *)malloc(sizeof(PointList));
      D1List = DataD1[i].list;
      D2List = DataD2[i].list;
    } else 
      DataD1[i].list = DataD2[i].list = NULL;
    DataD1[i].setName = 
      (char *)malloc(sizeof(char)*strlen(PlotData[i].setName));
    DataD2[i].setName = 
      (char *)malloc(sizeof(char)*strlen(PlotData[i].setName));
    strcpy(DataD1[i].setName,PlotData[i].setName);
    strcpy(DataD2[i].setName,PlotData[i].setName);
    while (theList) {
      int j;
 
      D1List->numPoints = D2List->numPoints = theList->numPoints;
      D1List->xvec = (double *)malloc(sizeof(double)*theList->numPoints);
      D1List->yvec = (double *)malloc(sizeof(double)*theList->numPoints);
      D1List->next = NULL;
      D2List->xvec = (double *)malloc(sizeof(double)*theList->numPoints);
      D2List->yvec = (double *)malloc(sizeof(double)*theList->numPoints);
      D2List->next = NULL;

      for (j=1;j<theList->numPoints-1;j++) {
        D1List->xvec[j] = D2List->xvec[j] = theList->xvec[j];
        D1List->yvec[j] = (theList->yvec[j+1] - theList->yvec[j-1]) /
                (theList->xvec[j+1] - theList->xvec[j-1]);
        D2List->yvec[j] = (theList->yvec[j+1] + theList->yvec[j-1] -
                 2.0*theList->yvec[j]) /
                ((theList->xvec[j+1] - theList->xvec[j]) * 
                 (theList->xvec[j]-theList->xvec[j-1]));
 
      }

      /* Extrapolate to get the endpoints ... */
      /* end near 0 */
      D1List->xvec[0] = theList->xvec[0];
      D1List->xvec[theList->numPoints-1] = 
       theList->xvec[theList->numPoints-1];
      m = (D1List->yvec[2]-D1List->yvec[1]) / 
          (theList->xvec[2]-theList->xvec[1]);
      b = D1List->yvec[1] - m*theList->xvec[1];
      D1List->yvec[0] = m*theList->xvec[0] + b;
      /* end near numPoints-1 */
      m = (D1List->yvec[theList->numPoints-2]-
           D1List->yvec[theList->numPoints-3]) / 
        (theList->xvec[theList->numPoints-2]-
         theList->xvec[theList->numPoints-3]);
      b = D1List->yvec [theList->numPoints-2] - 
          m*theList->xvec[theList->numPoints-2];
      D1List->yvec[theList->numPoints-1] = 
          m*theList->xvec[theList->numPoints-1] + b;

      /* Extrapolate to get the endpoints ... */
      /* end near 0 */
      D2List->xvec[0] = theList->xvec[0];
      D2List->xvec[theList->numPoints-1] = 
       theList->xvec[theList->numPoints-1];
      m = (D2List->yvec[2]-D2List->yvec[1]) / 
          (theList->xvec[2]-theList->xvec[1]);
      b = D2List->yvec[1] - m*theList->xvec[1];
      D2List->yvec[0] = m*theList->xvec[0] + b;
      /* end near numPoints-1 */
      m = (D2List->yvec[theList->numPoints-2]-
           D2List->yvec[theList->numPoints-3]) / 
        (theList->xvec[theList->numPoints-2]-
         theList->xvec[theList->numPoints-3]);
      b = D2List->yvec[theList->numPoints-2] - 
          m*theList->xvec[theList->numPoints-2];
      D2List->yvec[theList->numPoints-1] = 
           m*theList->xvec[theList->numPoints-1] + b;

      theList = theList->next;
      if (theList) {
        D1List->next = (PointList *)malloc(sizeof(PointList)); 
        D2List->next = (PointList *)malloc(sizeof(PointList)); 
        D1List = D1List->next;
        D2List = D2List->next;
      }
    }
  }
} 

void
Bounds(loX, loY, hiX, hiY, Ord)
  double *loX, *loY, *hiX, *hiY;
  int Ord;
{
  int i;
  PointList *theList;
  if ((Ord<1) || (Ord>2)) {
    printf ("Internal Error - Cannot eval deriv > 2 in Bounds.\n");
    exit(1);
  }
  *loX = *loY = *hiX = *hiY = 0.0;
  for (i=0;i<MAXSETS;i++) {
    if (Ord == 1)
      theList = DataD1[i].list;
    else
      theList = DataD2[i].list;
    while (theList) {
      int j;
 
      for (j=0;j<theList->numPoints;j++) {
        *loX = (theList->xvec[j]<*loX?theList->xvec[j]:*loX);
        *loY = (theList->yvec[j]<*loY?theList->yvec[j]:*loY);
        *hiX = (theList->xvec[j]>*hiX?theList->xvec[j]:*hiX);
        *hiY = (theList->yvec[j]>*hiY?theList->yvec[j]:*hiY);
      }
      theList = theList->next;
    }
  }
}