1 #include <math.h>
2 #include <R.h>
IntIndexSRC(double * L,double * R,double * p,double * q,int * N,int * M,int * Iindex,int * Mindex,int * Istrata,int * Mstrata)3 void IntIndexSRC(double *L,
4 double *R,
5 double *p,
6 double *q,
7 int *N,
8 int *M,
9 int *Iindex,
10 int *Mindex,
11 int *Istrata,
12 int *Mstrata){
13
14 int i,m,k,l;
15 k=0;
16
17 for (i=0; i<*N;i++){
18 for (m=0; m<*M;m++){
19 if ((L[i]==R[i] && p[m]==q[m] && L[i]==q[m]) /* point */
20 ||
21 (L[i]<q[m] && L[i]<=p[m] && R[i]>=q[m] && R[i]>p[m])) /* interval */
22 {
23 Iindex[k]=m+1;
24 k++;
25 }
26 }
27 Istrata[i]=k;
28 }
29 l=0;
30 for (m=0; m<*M;m++){
31 for (i=0; i<*N;i++){
32 if ((L[i]==R[i] && p[m]==q[m] && L[i]==q[m]) /* point */
33 ||
34 (L[i]<q[m] && L[i]<=p[m] && R[i]>=q[m] && R[i]>p[m])) /* interval */
35 {
36 Mindex[l]=i+1;
37 l++;
38 }
39 }
40 Mstrata[m]=l;
41 }
42 }
43
44
45
46
Turnb(int * Mstrata,int * Istrata,int * Mindex,int * Iindex,int * N,int * M,double * Z,double * nplme)47 void Turnb(int *Mstrata,
48 int *Istrata,
49 int *Mindex,
50 int *Iindex,
51 int *N,
52 int *M,
53 double *Z,
54 double *nplme){
55
56 int i,l,u,j,Iind, Mind;
57 double Ilast, ZI, ZM, Mlast, Zlast, ZMI;
58 Mlast=0;
59 for(i=0;i<*M;i++){
60
61 Zlast=0;
62 ZMI=0;
63
64 for(l=0;l<*N; l++){
65
66 Mlast=0;
67 ZM=0;
68 Mind=0;
69
70 for(u=Mstrata[l];u<Mstrata[l+1];u++){
71
72 Mind=Mindex[u];
73
74 Ilast=0;
75 ZI=0;
76 Iind=0;
77
78 for(j=Istrata[l]; j<Istrata[l+1];j++)
79 {
80 Iind=Iindex[j];
81 ZI=Z[Iind-1]+Ilast;
82 Ilast=ZI;
83 }
84
85 ZM=Z[Mind-1]/Ilast + Mlast;
86 Mlast=ZM;
87
88 }
89
90 ZMI=Mlast+Zlast;
91 Zlast=ZMI;
92 }
93
94 nplme[i]=Mlast;
95 }
96 }
97