← Back to branch summary

~ubuntu-branches/ubuntu/hardy/libterralib/hardy

« back to all changes in this revision

Viewing changes to src/terralib/application/TeSemivarModelFactory.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T Chen
  • Date: 2005-11-25 22:32:59 UTC
  • Revision ID: james.westby@ubuntu.com-20051125223259-3zubal8ux4ki4fjg
Tags: upstream-3.0.3b2
Import upstream version 3.0.3b2

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/terralib/application/TeSemivarModelFactory.cpp
 
1
/************************************************************************************
 
2
TerraLib - a library for developing GIS applications.
 
3
Copyright � 2001-2004 INPE and Tecgraf/PUC-Rio.
 
4
 
 
5
This code is part of the TerraLib library.
 
6
This library is free software; you can redistribute it and/or
 
7
modify it under the terms of the GNU Lesser General Public
 
8
License as published by the Free Software Foundation; either
 
9
version 2.1 of the License, or (at your option) any later version.
 
10
 
 
11
You should have received a copy of the GNU Lesser General Public
 
12
License along with this library.
 
13
 
 
14
The authors reassure the license terms regarding the warranties.
 
15
They specifically disclaim any warranties, including, but not limited to,
 
16
the implied warranties of merchantability and fitness for a particular purpose.
 
17
The library provided hereunder is on an "as is" basis, and the authors have no
 
18
obligation to provide maintenance, support, updates, enhancements, or modifications.
 
19
In no event shall INPE and Tecgraf / PUC-Rio be held liable to any party for direct,
 
20
indirect, special, incidental, or consequential damages arising out of the use
 
21
of this library and its documentation.
 
22
*************************************************************************************/
 
23
 
 
24
#include <TeSemivarModelFactory.h>
 
25
#include <math.h>
 
26
#include <vector>
 
27
using namespace std;
 
28
 
 
29
static TeEsfericSemivar_Factory mesf("Esferic");
 
30
static TeExponentialSemivar_Factory mexp("Exponential");
 
31
static TeGaussianSemivar_Factory mgau("Gaussian");
 
32
 
 
33
 
 
34
TeMatrix
 
35
TeEsfericSemivarModel :: calculate (TeMatrix& g)
 
36
{
 
37
        TeMatrix                mxesf;
 
38
        vector<double>  esf;
 
39
        int                             i, nlag;
 
40
        double                  c0, c1, a;
 
41
 
 
42
        c0=modeloefeitopepita_;
 
43
        c1=modelocontribuicao_;
 
44
        a=modeloalcance_;
 
45
        nlag = g.Nrow();
 
46
        mxesf.Init(nlag+1, 2);
 
47
 
 
48
        for (i=0; i<nlag; ++i)
 
49
        {
 
50
                if (g(i,0) == 0)
 
51
                        esf.push_back(c0);
 
52
 
 
53
                else if (g(i,0) > 0 && g(i,0) < 1000000.0)
 
54
                {
 
55
                        esf.push_back( c0 + c1*(1.5*(g(i,0)/a)) - 0.5*(pow(g(i,0)/a,3) ));
 
56
                        if (esf[i] > c0+c1)
 
57
                                esf[i]=c0+c1;
 
58
                }
 
59
 
 
60
                else 
 
61
            esf.push_back(c0+c1);
 
62
    
 
63
                // Escreve Saida
 
64
                mxesf(i+1,0) = g(i,0); // h
 
65
        mxesf(i+1,1) = esf[i]; // gama(h)
 
66
    }
 
67
 
 
68
        // quando h=0 => gama(h)=c0
 
69
        mxesf(0,1)=c0;
 
70
 
 
71
        // acrescentar o �ltimo valor
 
72
        mxesf(nlag,0)= mxesf(nlag-1,0) + (mxesf(nlag-1,0) - mxesf(nlag-2,0));
 
73
    mxesf(nlag,1)= mxesf(nlag-1,1);
 
74
 
 
75
        return mxesf;
 
76
}
 
77
 
 
78
TeMatrix
 
79
TeExponentialSemivarModel :: calculate (TeMatrix& g)
 
80
{
 
81
        TeMatrix                mxexp;
 
82
        vector<double>  expo;
 
83
        int                             i, nlag;
 
84
        double                  c0, c1, a;
 
85
 
 
86
        c0=modeloefeitopepita_;
 
87
        c1=modelocontribuicao_;
 
88
        a=modeloalcance_;
 
89
        nlag = g.Nrow();
 
90
        mxexp.Init(nlag+1, 2);
 
91
 
 
92
        for (i=0; i<nlag; ++i)
 
93
        {
 
94
                if (g(i,0) == 0)
 
95
                        expo.push_back(c0);
 
96
 
 
97
                //expo(i) = c0 + c1*( 1 - exp( -(3*g(i,1))/a) );
 
98
                expo.push_back(c0 + c1*( 1 - exp( -(3*g(i,0))/a) ));
 
99
 
 
100
                // Escreve Saida
 
101
                mxexp(i+1,0) = g(i,0); // h
 
102
        mxexp(i+1,1) = expo[i]; // gama(h)
 
103
    }
 
104
 
 
105
        // quando h=0 => gama(h)=c0
 
106
        mxexp(0,1)=c0;
 
107
 
 
108
        // acrescentar o �ltimo valor
 
109
        mxexp(nlag,0)= mxexp(nlag-1,0) + (mxexp(nlag-1,0) - mxexp(nlag-2,0));
 
110
    mxexp(nlag,1)= mxexp(nlag-1,1);
 
111
 
 
112
        return mxexp;
 
113
}
 
114
 
 
115
TeMatrix
 
116
TeGaussianSemivarModel :: calculate (TeMatrix& g)
 
117
{
 
118
        TeMatrix                mxgau;
 
119
        vector<double>  gau;
 
120
        int                             i, nlag;
 
121
        double                  c0, c1, a;
 
122
 
 
123
        c0=modeloefeitopepita_;
 
124
        c1=modelocontribuicao_;
 
125
        a=modeloalcance_;
 
126
        nlag = g.Nrow();
 
127
        mxgau.Init(nlag+1, 2);
 
128
 
 
129
        for (i=0; i<nlag; ++i)
 
130
        {
 
131
                if (g(i,0) == 0)
 
132
                        gau.push_back(c0);
 
133
 
 
134
           // gaus(i) = c0 + c1*( 1 - exp(-3*(g(i,1)/a)^2) );
 
135
                gau.push_back( c0 + c1*( 1 - exp(-3*pow(g(i,0)/a,2))) );
 
136
 
 
137
                // Escreve Saida
 
138
                mxgau(i+1,0) = g(i,0); // h
 
139
        mxgau(i+1,1) = gau[i]; // gama(h)
 
140
    }
 
141
 
 
142
        // quando h=0 => gama(h)=c0
 
143
        mxgau(0,1)=c0;
 
144
 
 
145
        // acrescentar o �ltimo valor
 
146
        mxgau(nlag,0)= mxgau(nlag-1,0) + (mxgau(nlag-1,0) - mxgau(nlag-2,0));
 
147
    mxgau(nlag,1)= mxgau(nlag-1,1);
 
148
 
 
149
        return mxgau;
 
150
}
 
151
 
 
152
 
 
153
/*
 
154
for i=1:nlag;
 
155
    if (g(i,1) == 0)
 
156
        sphe(i) = c0;
 
157
        gaus(i) = c0;
 
158
        expo(i) = c0;
 
159
    elseif (g(i,1) > 0 & g(i,1) < a)
 
160
        sphe(i) = c0 + c1*((1.5*(g(i,1)/a) - 0.5*((g(i,1)/a)^3)));
 
161
    else 
 
162
        if (g(i,1) >= a)
 
163
            sphe(i)=c0+c1;
 
164
        end
 
165
    end
 
166
    gaus(i) = c0 + c1*( 1 - exp(-3*(g(i,1)/a)^2) );
 
167
    expo(i) = c0 + c1*( 1 - exp( -(3*g(i,1))/a) );
 
168
    
 
169
    % Escreve Saida
 
170
    m(i+1,1) = g(i,1);% lag distance
 
171
 
 
172
    switch (model)
 
173
        case 1
 
174
            m(i+1,2) = sphe(i);
 
175
        case 2
 
176
            m(i+1,2) = expo(i);
 
177
        case 3
 
178
            m(i+1,2) = gaus(i);
 
179
        otherwise      
 
180
            m(i+1,2) = 0;
 
181
    end
 
182
end
 
183
 
 
184
%Para o modelo quando h=0 => gamma(h)=c0
 
185
m(1,2)=c0;
 
186
 
 
187
%Para o modelo quando h >= a => gamma(h)=c0+c1
 
188
switch (model)
 
189
    case 1 %Para modelo esferico
 
190
    m(nlag+2,2)= c0+c1;
 
191
 
 
192
    case 2 %Para modelo exponencial
 
193
    m(nlag+2,2)= m(nlag+1,2);
 
194
    
 
195
    case 3 %Para modelo gaussiano
 
196
    m(nlag+2,2)= m(nlag+1,2);
 
197
end
 
198
*/