00001
00002
00003
00004
00005
00006
00007
00008
00009 #include <stdlib.h>
00010 #include "tablefunct.h"
00011 #include "intools.h"
00012
00013
00014
00015
00016
00017
00018
00019
00020 tablefunct::tablefunct ()
00021 {
00022 itype = piecewiselin;
00023 asize=0;
00024 x=NULL;
00025 y=NULL;
00026 file1 = new char[1000];
00027 }
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038 tablefunct::tablefunct (long n)
00039 {
00040 itype = piecewiselin;
00041 asize = n;
00042 x = new double [n];
00043 y = new double [n];
00044 file1 = new char[1000];
00045 }
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056 tablefunct::tablefunct (long n,long it)
00057 {
00058 switch (it){
00059 case 1:{
00060 itype = piecewiselin;
00061 break;
00062 }
00063 case 2:{
00064 itype = piecewiseconst;
00065 break;
00066 }
00067 case 3:{
00068 itype = lagrange;
00069 break;
00070 }
00071 case 4:{
00072 itype = piecewiselin2;
00073 break;
00074 }
00075 case 5:{
00076 itype = dirac2;
00077 break;
00078 }
00079 default:{
00080 print_err ("unknown type of interpolation is required",__FILE__,__LINE__,__func__);
00081 }
00082 }
00083
00084 asize = n;
00085 x = new double [n];
00086 y = new double [n];
00087 }
00088
00089
00090
00091
00092
00093
00094
00095
00096 tablefunct::~tablefunct()
00097 {
00098 delete [] x;
00099 delete [] y;
00100 delete [] file1;
00101 }
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114 void tablefunct::read(XFILE *in)
00115 {
00116
00117
00118
00119
00120
00121 xfscanf (in,"%k%m%k%ld","approx_type", &interpoltype_kwdset, (int*)&itype, "ntab_items", &asize);
00122
00123 x=new double[asize];
00124 y=new double[asize];
00125
00126 readval (in);
00127 }
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138 void tablefunct::read_data_file(XFILE *in)
00139 {
00140 XFILE *in_data;
00141
00142 xfscanf(in, " %a", file1);
00143
00144 in_data = xfopen(file1,"r");
00145
00146
00147
00148
00149
00150 xfscanf (in_data,"%k%m%k%ld","approx_type", &interpoltype_kwdset, (int*)&itype, "ntab_items", &asize);
00151
00152 x=new double[asize];
00153 y=new double[asize];
00154
00155 readval (in_data);
00156
00157 xfclose(in_data);
00158 }
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168 void tablefunct::readval (XFILE *in)
00169 {
00170 long i;
00171 for (i=0;i<asize;i++){
00172 xfscanf (in,"%lf",&x[i]);
00173 xfscanf (in,"%lf",&y[i]);
00174 }
00175
00176
00177
00178
00179 }
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191 void tablefunct::print(FILE *out)
00192 {
00193 long i;
00194 fprintf (out,"\n %d %ld",itype,asize);
00195 for (i=0;i<asize;i++){
00196 fprintf (out,"\n %le ",x[i]);
00197 fprintf (out,"%le ",y[i]);
00198 }
00199 }
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209 void tablefunct::print_data_file(FILE *out)
00210 {
00211 fprintf (out,"\n %s",file1);
00212 }
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224 void tablefunct::read_prop(FILE *in)
00225 {
00226 long i;
00227 getint(in, (int&)itype);
00228 getlong(in, asize);
00229 x=new double[asize];
00230 y=new double[asize];
00231 for (i=0;i<asize;i++)
00232 {
00233 getdouble (in, x[i]);
00234 getdouble (in, y[i]);
00235 }
00236 }
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248 void tablefunct::datacheck ()
00249 {
00250 long i;
00251
00252 for (i=0;i<asize-1;i++){
00253 if (x[i]<x[i+1]){
00254 continue;
00255 }
00256 else{
00257 print_err("the %ld-th component (%lf, %lf) does not satisfy condition x[0] < x[1] < x[2] < ... <x[asize-1]", __FILE__, __LINE__, __func__,i,x[i],x[i+1]);
00258 }
00259 }
00260 }
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278 double tablefunct::piecewise_const_interpol (double temp)
00279 {
00280 long i;
00281 double f;
00282
00283 if (temp<x[0] || temp>=x[asize-1]){
00284 print_err("required value %le is out of range <%le;%le>", __FILE__, __LINE__, __func__, temp, x[0], x[asize-1]);
00285 exit(2);
00286 }
00287 for (i=0;i<asize-1;i++){
00288 if (x[i]<=temp && temp<x[i+1]){
00289 f=y[i];
00290 break;
00291 }
00292 }
00293 return f;
00294 }
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306 double tablefunct::diracinterpol2(double temp)
00307 {
00308 long i;
00309 double iv;
00310 for (i=0;i<asize;i++){
00311 if (i==asize-1){
00312 iv=y[i];
00313 break;
00314 }
00315 if ((i==0) && (x[i] >= temp))
00316 {
00317 iv=y[i];
00318 break;
00319 }
00320 if (temp>x[i] && temp<=x[i+1]){
00321 iv=y[i+1];
00322 break;
00323 }
00324 }
00325 return iv;
00326 }
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339 double tablefunct::diracinterpol2(double temp, double &k)
00340 {
00341 long i;
00342 double iv;
00343 for (i=0;i<asize;i++){
00344 if (i==asize-1){
00345 iv=y[i];
00346 break;
00347 }
00348 if ((i==0) && (x[i] >= temp))
00349 {
00350 iv=y[i];
00351 break;
00352 }
00353 if (temp>x[i] && temp<=x[i+1]){
00354 iv=y[i+1];
00355 break;
00356 }
00357 }
00358 k = 0.0;
00359 return iv;
00360 }
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375 double tablefunct::piecewise_linear_interpol (double temp)
00376 {
00377 long i;
00378 double f,k;
00379
00380 if (temp<x[0] || temp>x[asize-1]){
00381 print_err("required value %le is out of range <%le;%le>", __FILE__, __LINE__, __func__, temp, x[0], x[asize-1]);
00382 exit(2);
00383 }
00384
00385 for (i=0;i<asize-1;i++){
00386 if (x[i]<=temp && temp<=x[i+1]){
00387 k =(y[i+1]-y[i])/(x[i+1]-x[i]);
00388 f=y[i]+(temp-x[i])*k;
00389 break;
00390 }
00391 }
00392
00393 return f;
00394 }
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408 double tablefunct::inverse_piecewise_linear_interpol (double temp)
00409 {
00410 long i;
00411 double f,k;
00412
00413 if (temp<y[0] || temp>y[asize-1]){
00414 print_err("required value %le is out of range <%le;%le>", __FILE__, __LINE__, __func__, temp, y[0], y[asize-1]);
00415 exit(2);
00416 }
00417
00418 for (i=0;i<asize-1;i++){
00419 if (y[i]<=temp && temp<=y[i+1]){
00420 k =(x[i+1]-x[i])/(y[i+1]-y[i]);
00421 f=x[i]+(temp-y[i])*k;
00422 break;
00423 }
00424 }
00425
00426 return f;
00427 }
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438 double tablefunct::lininterpol2 (double temp, double &k)
00439 {
00440
00441
00442 double q,iv;
00443
00444 int i,h=0;
00445
00446 for (i=0;i<asize;i++){
00447 if ((x[i] >= temp) || (i == asize-1)){
00448 if (i == 0)
00449 {
00450 k =(y[i+1]-y[i])/(x[i+1]-x[i]);
00451 q =y[i]-k*x[i];
00452 iv=k*temp+q;
00453 h=1;
00454 break;
00455 }
00456 k =(y[i]-y[i-1])/(x[i]-x[i-1]);
00457 q =y[i]-k*x[i];
00458 iv=k*temp+q;
00459 h=1;
00460 break;
00461 }
00462 }
00463 if(h==0){
00464 print_err("out of range value is required", __FILE__, __LINE__, __func__);
00465 exit(2);
00466 }
00467
00468 return iv;
00469 }
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482 double tablefunct::lininterpol3 (double temp, double &k)
00483 {
00484
00485
00486 double q,iv;
00487
00488 int i,h=0;
00489
00490 for (i=0;i<asize;i++){
00491 if ((y[i] >= temp) || (i == asize-1)){
00492 if (i == 0)
00493 {
00494 k =(x[i+1]-x[i])/(y[i+1]-y[i]);
00495 q =x[i]-k*y[i];
00496 iv=k*temp+q;
00497 h=1;
00498 break;
00499 }
00500 k =(x[i]-x[i-1])/(y[i]-y[i-1]);
00501 q =x[i]-k*y[i];
00502 iv=k*temp+q;
00503 h=1;
00504 break;
00505 }
00506 }
00507 if(h==0){
00508 print_err("out of range value is required", __FILE__, __LINE__, __func__);
00509 exit(2);
00510 }
00511
00512 return iv;
00513 }
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524 double tablefunct::laginterpol(double temp)
00525 {
00526 double y1=0.,nom,denom;
00527 long i,j;
00528
00529 for(i=0; i<asize; i++){
00530 nom=y[i];
00531 denom=1.;
00532 for(j=0; j<asize; j++){
00533 if(i==j) continue;
00534 nom=nom*(temp-x[j]);
00535 denom=denom*(x[i]-x[j]);
00536 if(denom<1.0e-8) continue;
00537 }
00538 y1=y1+nom/denom;
00539 }
00540
00541 return y1;
00542 }
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553 double tablefunct::derivative (double temp)
00554 {
00555 long i;
00556 double d;
00557
00558 if (temp<x[0] || temp>x[asize-1]){
00559 print_err("required value %le is out of range <%le, %le>", __FILE__, __LINE__, __func__, temp, x[0], x[asize-1]);
00560 exit(2);
00561 }
00562
00563 for (i=0;i<asize-1;i++){
00564 if (x[i]<=temp && temp<=x[i+1]){
00565 d=(y[i+1]-y[i])/(x[i+1]-x[i]);
00566 break;
00567 }
00568 }
00569
00570 return d;
00571 }
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581 double tablefunct::inv_derivative (double temp)
00582 {
00583 long i;
00584 double d;
00585
00586 if (temp<y[0] || temp>y[asize-1]){
00587 print_err("required value is out of range", __FILE__, __LINE__, __func__);
00588 exit(2);
00589 }
00590
00591 for (i=0;i<asize-1;i++){
00592 if (y[i]<=temp && temp<=y[i+1]){
00593 d=(x[i+1]-x[i])/(y[i+1]-y[i]);
00594 break;
00595 }
00596 }
00597
00598 return d;
00599 }
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611 double tablefunct::getval(double temp)
00612 {
00613
00614 double iv;
00615 if (asize<=1 || x==NULL || y==NULL){
00616 print_err("wrong size of arrays or wrong number of components", __FILE__, __LINE__, __func__);
00617 exit(2);
00618 }
00619
00620 switch(itype){
00621 case piecewiselin:{
00622 iv=piecewise_linear_interpol (temp);
00623 break;
00624 }
00625 case piecewiseconst:{
00626 iv=piecewise_const_interpol (temp);
00627 break;
00628 }
00629 case dirac2:{
00630 double dummy;
00631 iv=diracinterpol2(temp, dummy);
00632 break;
00633 }
00634 case lagrange:{
00635 iv=laginterpol(temp);
00636 break;
00637 }
00638 case piecewiselin2:{
00639 iv=piecewise_linear_interpol (temp);
00640 break;
00641 }
00642 default:{
00643 print_err("unknown type of interpolation", __FILE__, __LINE__, __func__);
00644 exit(2);
00645 }
00646 }
00647
00648 return iv;
00649 }
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660 double tablefunct::getinvval(double temp)
00661 {
00662
00663 double iv;
00664 if (asize<=1 || x==NULL || y==NULL){
00665 print_err("wrong size of arrays or wrong number of components", __FILE__, __LINE__, __func__);
00666 exit(2);
00667 }
00668
00669 switch(itype){
00670 case piecewiselin:{
00671 iv=inverse_piecewise_linear_interpol (temp);
00672 break;
00673 }
00674 default:{
00675 print_err("unknown type of interpolation", __FILE__, __LINE__, __func__);
00676 exit(2);
00677 }
00678 }
00679
00680 return iv;
00681 }
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693 double tablefunct::getval2(double temp, double &k)
00694 {
00695
00696 double iv;
00697 if (asize<=1 || x==NULL || y==NULL){
00698 fprintf (stderr,"\ntablefunct::getval2\nWrong size of the array (%s, line %d)",__FILE__,__LINE__);
00699 exit(2);
00700 }
00701
00702 switch(itype){
00703 case piecewiselin2:{
00704 iv=lininterpol2(temp, k);
00705 break;
00706 }
00707 case dirac2:{
00708 iv=diracinterpol2(temp, k);
00709 break;
00710 }
00711 default:{
00712 fprintf (stderr,"\ntablefunct::getval2\nUnknown type of interpolation (%s, line %d)",__FILE__,__LINE__);
00713 exit(2);
00714 }
00715 }
00716
00717 return iv;
00718 }
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731 double tablefunct::getval3(double temp, double &k)
00732 {
00733
00734 double iv;
00735 if (asize<=1 || x==NULL || y==NULL){
00736 print_err("unknown definition of material parameter is required",__FILE__,__LINE__,__func__);
00737 exit(2);
00738 }
00739
00740 switch(itype){
00741 case piecewiselin:
00742 case piecewiselin2:{
00743 iv=lininterpol3(temp, k);
00744 break;
00745 }
00746 default:{
00747 print_err("unknown definition of material parameter is required",__FILE__,__LINE__,__func__);
00748 exit(2);
00749 }
00750 }
00751
00752 return iv;
00753 }
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764 double tablefunct::getderiv (double temp)
00765 {
00766 double d;
00767
00768 if (asize<=1 || x==NULL || y==NULL){
00769 print_err("wrong size of arrays or wrong number of components", __FILE__, __LINE__, __func__);
00770 exit(2);
00771 }
00772
00773 switch(itype){
00774 case piecewiselin:{
00775 d = derivative (temp);
00776 break;
00777 }
00778 default:{
00779 print_err("unknown type of interpolation", __FILE__, __LINE__, __func__);
00780 exit(2);
00781 }
00782 }
00783
00784 return d;
00785 }
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798 void tablefunct::copy (tablefunct &tf)
00799 {
00800 long i;
00801
00802
00803 itype = tf.itype;
00804
00805 asize = tf.asize;
00806
00807 if (x!=NULL)
00808 delete [] x;
00809 x = new double [asize];
00810
00811 if (y!=NULL)
00812 delete [] y;
00813 y = new double [asize];
00814
00815 for (i=0;i<asize;i++){
00816 x[i]=tf.x[i];
00817 y[i]=tf.y[i];
00818 }
00819 }
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833 long tablefunct::compare(tablefunct &tf)
00834 {
00835 long i;
00836
00837
00838 if (itype != tf.itype)
00839 return 1;
00840
00841
00842 if (asize != tf.asize)
00843 return 1;
00844
00845
00846 for (i=0;i<asize;i++){
00847 if (x[i] != tf.x[i])
00848 return 1;
00849 if (y[i] != tf.y[i])
00850 return 1;
00851 }
00852 return 0;
00853 }