1 // GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege
2 //
3 // See the LICENSE.txt file for license information. Please report all
4 // issues on https://gitlab.onelab.info/getdp/getdp/issues.
5 //
6 // Contributor(s):
7 // Johan Gyselinck
8 // Ruth Sabariego
9
10 #include <math.h>
11 #include "F.h"
12 #include "GeoData.h"
13 #include "DofData.h"
14 #include "Cal_Value.h"
15 #include "Message.h"
16
17 extern struct CurrentData Current ;
18
19 #define SQU(a) ((a)*(a))
20 #define TWO_PI 6.2831853071795865
21
22 /* ------------------------------------------------------------------------ */
23 /* Simple Extended Math */
24 /* ------------------------------------------------------------------------ */
25
F_Hypot(F_ARG)26 void F_Hypot(F_ARG)
27 {
28 int k;
29 double tmp;
30
31 if(A->Type != SCALAR || (A+1)->Type != SCALAR)
32 Message::Error("Non scalar argument(s) for function 'Hypot'");
33
34 if (Current.NbrHar == 1){
35 V->Val[0] = sqrt(A->Val[0]*A->Val[0]+(A+1)->Val[0]*(A+1)->Val[0]) ;
36 }
37 else {
38 tmp = sqrt(A->Val[0]*A->Val[0]+(A+1)->Val[0]*(A+1)->Val[0]) ;
39 for (k = 0 ; k < Current.NbrHar ; k += 2) {
40 V->Val[MAX_DIM* k ] = tmp ;
41 V->Val[MAX_DIM*(k+1)] = 0. ;
42 }
43 }
44 V->Type = SCALAR;
45 }
46
F_TanhC2(F_ARG)47 void F_TanhC2(F_ARG)
48 {
49 // @Jon: check this and the behavior of the overall function for large args
50 // lim_x\to\inf cosh(x) = +\inf
51 // lim_x\to\inf sinh(x) = 1
52
53 double denom =
54 SQU(cosh(A->Val[0])*cos(A->Val[MAX_DIM])) +
55 SQU(sinh(A->Val[0])*sin(A->Val[MAX_DIM]));
56 //printf("arg=%g cosh(arg)=%g\n", A->Val[0], cosh(A->Val[0]));
57
58
59 V->Val[0] = sinh(A->Val[0])*cosh(A->Val[0]) / denom ;
60 V->Val[MAX_DIM] = sin(A->Val[MAX_DIM])*cos(A->Val[MAX_DIM]) / denom ;
61 V->Type = SCALAR ;
62
63 /* printf("numer_real=%g, numer_imag=%g, denom = %g\n",
64 sinh(A->Val[0])*cosh(A->Val[0]) ,
65 sin(A->Val[MAX_DIM])*cos(A->Val[MAX_DIM]),
66 denom); */
67
68 }
69
70 /* ------------------------------------------------------------------------ */
71 /* General Tensor Functions */
72 /* ------------------------------------------------------------------------ */
73
F_Transpose(F_ARG)74 void F_Transpose(F_ARG)
75 {
76 if(A->Type != TENSOR_DIAG && A->Type != TENSOR_SYM && A->Type != TENSOR)
77 Message::Error("Wrong type of argument for function 'Transpose'");
78
79 Cal_TransposeValue(A,V);
80 }
81
F_Inv(F_ARG)82 void F_Inv(F_ARG)
83 {
84 if(A->Type != TENSOR_DIAG && A->Type != TENSOR_SYM && A->Type != TENSOR)
85 Message::Error("Wrong type of argument for function 'Inverse'");
86
87 Cal_InvertValue(A,V);
88 }
89
F_Det(F_ARG)90 void F_Det(F_ARG)
91 {
92 if(A->Type != TENSOR_DIAG && A->Type != TENSOR_SYM && A->Type != TENSOR)
93 Message::Error("Wrong type of argument for function 'Det'");
94
95 Cal_DetValue(A,V);
96 }
97
F_Trace(F_ARG)98 void F_Trace(F_ARG)
99 {
100 if(A->Type != TENSOR_DIAG && A->Type != TENSOR_SYM && A->Type != TENSOR)
101 Message::Error("Wrong type of argument for function 'Trace'");
102
103 Cal_TraceValue(A,V);
104 }
105
F_RotateXYZ(F_ARG)106 void F_RotateXYZ(F_ARG)
107 {
108 // Apply a (X_1 Y_2 Z_3) rotation matrix using Euler (Tait-Bryan) angles
109 double ca, sa, cb, sb, cc, sc ;
110 struct Value Rot ;
111
112 if((A->Type != TENSOR_DIAG && A->Type != TENSOR_SYM && A->Type != TENSOR &&
113 A->Type != VECTOR) ||
114 (A+1)->Type != SCALAR || (A+2)->Type != SCALAR || (A+3)->Type != SCALAR)
115 Message::Error("Wrong type of argument(s) for function 'Rotate'");
116
117 ca = cos((A+1)->Val[0]) ; sa = sin((A+1)->Val[0]) ;
118 cb = cos((A+2)->Val[0]) ; sb = sin((A+2)->Val[0]) ;
119 cc = cos((A+3)->Val[0]) ; sc = sin((A+3)->Val[0]) ;
120
121 Rot.Type = TENSOR ;
122 Cal_ZeroValue(&Rot);
123 Rot.Val[0] = cb*cc; Rot.Val[1] = -cb*sc; Rot.Val[2] = sb;
124 Rot.Val[3] = sa*sb*cc+ca*sc; Rot.Val[4] = -sa*sb*sc+ca*cc; Rot.Val[5] = -sa*cb;
125 Rot.Val[6] = -ca*sb*cc+sa*sc; Rot.Val[7] = ca*sb*sc+sa*cc; Rot.Val[8] = ca*cb;
126
127 Cal_RotateValue(&Rot,A,V);
128 }
129
130 /* ------------------------------------------------------------------------ */
131 /* Norm */
132 /* ------------------------------------------------------------------------ */
133
F_Norm(F_ARG)134 void F_Norm(F_ARG)
135 {
136 int k ;
137
138 switch(A->Type) {
139
140 case SCALAR :
141 if (Current.NbrHar == 1) {
142 V->Val[0] = fabs(A->Val[0]) ;
143 }
144 else {
145 for (k = 0 ; k < Current.NbrHar ; k += 2 ) {
146 V->Val[MAX_DIM*k] = sqrt(SQU(A->Val[MAX_DIM*k]) +
147 SQU(A->Val[MAX_DIM*(k+1)]));
148 V->Val[MAX_DIM*(k+1)] = 0. ;
149 }
150 }
151 break ;
152
153 case VECTOR :
154 if (Current.NbrHar == 1) {
155 V->Val[0] = sqrt(SQU(A->Val[0]) + SQU(A->Val[1]) + SQU(A->Val[2])) ;
156 }
157 else {
158 for (k = 0 ; k < Current.NbrHar ; k += 2 ) {
159 V->Val[MAX_DIM*k] = sqrt(SQU(A->Val[MAX_DIM* k ]) +
160 SQU(A->Val[MAX_DIM* k +1]) +
161 SQU(A->Val[MAX_DIM* k +2]) +
162 SQU(A->Val[MAX_DIM*(k+1) ]) +
163 SQU(A->Val[MAX_DIM*(k+1)+1]) +
164 SQU(A->Val[MAX_DIM*(k+1)+2])) ;
165 V->Val[MAX_DIM*(k+1)] = 0. ;
166 }
167 }
168 break ;
169
170 default :
171 Message::Error("Wrong type of argument for function 'Norm'");
172 break;
173 }
174
175 V->Type = SCALAR ;
176 }
177
178 /* ------------------------------------------------------------------------ */
179 /* Square Norm */
180 /* ------------------------------------------------------------------------ */
181
F_SquNorm(F_ARG)182 void F_SquNorm(F_ARG)
183 {
184 int k ;
185
186 switch(A->Type) {
187
188 case SCALAR :
189 if (Current.NbrHar == 1) {
190 V->Val[0] = SQU(A->Val[0]) ;
191 }
192 else {
193 for (k = 0 ; k < Current.NbrHar ; k += 2 ) {
194 V->Val[MAX_DIM*k] = SQU(A->Val[MAX_DIM*k]) + SQU(A->Val[MAX_DIM*(k+1)]);
195 V->Val[MAX_DIM*(k+1)] = 0. ;
196 }
197 }
198 break ;
199
200 case VECTOR :
201 if (Current.NbrHar == 1) {
202 V->Val[0] = SQU(A->Val[0]) + SQU(A->Val[1]) + SQU(A->Val[2]) ;
203 }
204 else {
205 for (k = 0 ; k < Current.NbrHar ; k += 2 ) {
206 V->Val[MAX_DIM*k] = SQU(A->Val[MAX_DIM* k ]) +
207 SQU(A->Val[MAX_DIM* k +1]) +
208 SQU(A->Val[MAX_DIM* k +2]) +
209 SQU(A->Val[MAX_DIM*(k+1) ]) +
210 SQU(A->Val[MAX_DIM*(k+1)+1]) +
211 SQU(A->Val[MAX_DIM*(k+1)+2]) ;
212 V->Val[MAX_DIM*(k+1)] = 0. ;
213 }
214 }
215 break ;
216
217 default :
218 Message::Error("Wrong type of argument for function 'SquNorm'");
219 break;
220 }
221
222 V->Type = SCALAR ;
223 }
224
225 /* ------------------------------------------------------------------------ */
226 /* Unit */
227 /* ------------------------------------------------------------------------ */
228
F_Unit(F_ARG)229 void F_Unit(F_ARG)
230 {
231 int k ;
232 double Norm ;
233
234 switch(A->Type) {
235
236 case SCALAR :
237 if (Current.NbrHar == 1) {
238 V->Val[0] = 1. ;
239 }
240 else {
241 for (k = 0 ; k < Current.NbrHar ; k += 2 ) {
242 V->Val[MAX_DIM* k ] = 1. ;
243 V->Val[MAX_DIM*(k+1)] = 0. ;
244 }
245 }
246 V->Type = SCALAR ;
247 break ;
248
249 case VECTOR :
250 if (Current.NbrHar == 1) {
251 Norm = sqrt(SQU(A->Val[0]) + SQU(A->Val[1]) + SQU(A->Val[2])) ;
252 if (Norm > 1.e-30) { /* Attention: tolerance */
253 V->Val[0] = A->Val[0]/Norm ;
254 V->Val[1] = A->Val[1]/Norm ;
255 V->Val[2] = A->Val[2]/Norm ;
256 }
257 else {
258 V->Val[0] = 0. ;
259 V->Val[1] = 0. ;
260 V->Val[2] = 0. ;
261 }
262 }
263 else {
264 for (k = 0 ; k < Current.NbrHar ; k += 2 ) {
265 Norm = sqrt(SQU(A->Val[MAX_DIM* k ]) +
266 SQU(A->Val[MAX_DIM* k +1]) +
267 SQU(A->Val[MAX_DIM* k +2]) +
268 SQU(A->Val[MAX_DIM*(k+1) ]) +
269 SQU(A->Val[MAX_DIM*(k+1)+1]) +
270 SQU(A->Val[MAX_DIM*(k+1)+2])) ;
271 if (Norm > 1.e-30) { /* Attention: tolerance */
272 V->Val[MAX_DIM* k ] = A->Val[MAX_DIM* k ]/Norm ;
273 V->Val[MAX_DIM* k +1] = A->Val[MAX_DIM* k +1]/Norm ;
274 V->Val[MAX_DIM* k +2] = A->Val[MAX_DIM* k +2]/Norm ;
275 V->Val[MAX_DIM*(k+1) ] = A->Val[MAX_DIM*(k+1) ]/Norm ;
276 V->Val[MAX_DIM*(k+1)+1] = A->Val[MAX_DIM*(k+1)+1]/Norm ;
277 V->Val[MAX_DIM*(k+1)+2] = A->Val[MAX_DIM*(k+1)+2]/Norm ;
278 }
279 else {
280 V->Val[MAX_DIM* k ] = 0 ;
281 V->Val[MAX_DIM* k +1] = 0 ;
282 V->Val[MAX_DIM* k +2] = 0 ;
283 V->Val[MAX_DIM*(k+1) ] = 0 ;
284 V->Val[MAX_DIM*(k+1)+1] = 0 ;
285 V->Val[MAX_DIM*(k+1)+2] = 0 ;
286 }
287 }
288 }
289 V->Type = VECTOR ;
290 break ;
291
292 default :
293 Message::Error("Wrong type of argument for function 'Unit'");
294 break;
295 }
296 }
297
298 /* ------------------------------------------------------------------------ */
299 /* ScalarUnit */
300 /* ------------------------------------------------------------------------ */
301
F_ScalarUnit(F_ARG)302 void F_ScalarUnit(F_ARG)
303 {
304 int k ;
305
306 if (Current.NbrHar == 1) {
307 V->Val[0] = 1. ;
308 }
309 else {
310 for (k = 0 ; k < Current.NbrHar ; k += 2 ) {
311 V->Val[MAX_DIM* k ] = 1. ;
312 V->Val[MAX_DIM*(k+1)] = 0. ;
313 }
314 }
315 V->Type = SCALAR ;
316 }
317
318 /* ------------------------------------------------------------------------ */
319 /* Time Functions */
320 /* ------------------------------------------------------------------------ */
321
322 /* Interesting only because it allows the same formal expression in both
323 Time and Frequency domains ! */
324
325 /* cos ( w * $Time + phi ) */
326
F_Cos_wt_p(F_ARG)327 void F_Cos_wt_p (F_ARG)
328 {
329 if (Current.NbrHar == 1)
330 V->Val[0] = cos(Fct->Para[0] * Current.Time + Fct->Para[1]) ;
331 else if (Current.NbrHar == 2) {
332 V->Val[0] = cos(Fct->Para[1]) ;
333 V->Val[MAX_DIM] = sin(Fct->Para[1]) ;
334 }
335 else {
336 Message::Error("Too many harmonics for function 'Cos_wt_p'") ;
337 }
338 V->Type = SCALAR ;
339 }
340
341 /* sin ( w * $Time + phi ) */
342
F_Sin_wt_p(F_ARG)343 void F_Sin_wt_p (F_ARG)
344 {
345 if (Current.NbrHar == 1)
346 V->Val[0] = sin(Fct->Para[0] * Current.Time + Fct->Para[1]) ;
347 else if (Current.NbrHar == 2){
348 V->Val[0] = sin(Fct->Para[1]) ;
349 V->Val[MAX_DIM] = -cos(Fct->Para[1]) ;
350 }
351 else {
352 Message::Error("Too many harmonics for function 'Sin_wt_p'") ;
353 }
354 V->Type = SCALAR ;
355 }
356
F_Complex_MH(F_ARG)357 void F_Complex_MH(F_ARG)
358 {
359 int NbrFreq, NbrComp, i, j, k, l ;
360 struct Value R;
361 double * Val_Pulsation ;
362
363 NbrFreq = Fct->NbrParameters ;
364 NbrComp = Fct->NbrArguments ;
365 if (NbrComp != 2*NbrFreq)
366 Message::Error("Number of components does not equal twice the number "
367 "of frequencies in Complex_MH") ;
368
369 R.Type = A->Type ;
370 Cal_ZeroValue(&R);
371
372 if (Current.NbrHar != 1) {
373 Val_Pulsation = Current.DofData->Val_Pulsation ;
374 for (i=0 ; i<NbrFreq ; i++) {
375 for (j = 0 ; j < Current.NbrHar/2 ; j++)
376 if (fabs (Val_Pulsation[j]-TWO_PI*Fct->Para[i]) <= 1e-10 * Val_Pulsation[j]) {
377 for (k=2*j,l=2*i ; k<2*j+2 ; k++,l++) {
378 switch(A->Type){
379 case SCALAR :
380 R.Val[MAX_DIM*k ] += (A+l)->Val[0] ;
381 break;
382 case VECTOR :
383 case TENSOR_DIAG :
384 R.Val[MAX_DIM*k ] += (A+l)->Val[0] ;
385 R.Val[MAX_DIM*k+1] += (A+l)->Val[1] ;
386 R.Val[MAX_DIM*k+2] += (A+l)->Val[2] ;
387 break;
388 case TENSOR_SYM :
389 R.Val[MAX_DIM*k ] += (A+l)->Val[0] ;
390 R.Val[MAX_DIM*k+1] += (A+l)->Val[1] ;
391 R.Val[MAX_DIM*k+2] += (A+l)->Val[2] ;
392 R.Val[MAX_DIM*k+3] += (A+l)->Val[3] ;
393 R.Val[MAX_DIM*k+4] += (A+l)->Val[4] ;
394 R.Val[MAX_DIM*k+5] += (A+l)->Val[5] ;
395 break;
396 case TENSOR :
397 R.Val[MAX_DIM*k ] += (A+l)->Val[0] ;
398 R.Val[MAX_DIM*k+1] += (A+l)->Val[1] ;
399 R.Val[MAX_DIM*k+2] += (A+l)->Val[2] ;
400 R.Val[MAX_DIM*k+3] += (A+l)->Val[3] ;
401 R.Val[MAX_DIM*k+4] += (A+l)->Val[4] ;
402 R.Val[MAX_DIM*k+5] += (A+l)->Val[5] ;
403 R.Val[MAX_DIM*k+6] += (A+l)->Val[6] ;
404 R.Val[MAX_DIM*k+7] += (A+l)->Val[7] ;
405 R.Val[MAX_DIM*k+8] += (A+l)->Val[8] ;
406 break;
407 default :
408 Message::Error("Unknown type of arguments in function 'Complex_MH'");
409 break;
410 }
411 }
412 }
413 }
414 } else { /* time domain */
415 for (i=0 ; i<NbrFreq ; i++) {
416 Cal_AddMultValue (&R, A+2*i, cos(TWO_PI*Fct->Para[i]*Current.Time), &R) ;
417 Cal_AddMultValue (&R, A+2*i+1, -sin(TWO_PI*Fct->Para[i]*Current.Time), &R) ;
418 }
419 }
420 Cal_CopyValue(&R,V);
421 }
422
423 /* ------------------------------------------------------------------------ */
424 /* Period */
425 /* ------------------------------------------------------------------------ */
426
F_Period(F_ARG)427 void F_Period (F_ARG)
428 {
429 if (Current.NbrHar == 1)
430 V->Val[0] = fmod(A->Val[0], Fct->Para[0])
431 + ((A->Val[0] < 0.)? Fct->Para[0] : 0.) ;
432 else
433 Message::Error("Function 'F_Period' not valid for Complex");
434 V->Type = SCALAR ;
435 }
436
437 /* ------------------------------------------------------------------------ */
438 /* Interval */
439 /* ------------------------------------------------------------------------ */
440
F_Interval(F_ARG)441 void F_Interval (F_ARG)
442 {
443 int k;
444 double tmp;
445
446 if (Current.NbrHar == 1) {
447 V->Val[0] =
448 A->Val[0] > (A+1)->Val[0] + Fct->Para[0] * Fct->Para[2] &&
449 A->Val[0] < (A+2)->Val[0] + Fct->Para[1] * Fct->Para[2] ;
450 }
451 else {
452 tmp =
453 A->Val[0] > (A+1)->Val[0] + Fct->Para[0] * Fct->Para[2] &&
454 A->Val[0] < (A+2)->Val[0] + Fct->Para[1] * Fct->Para[2] ;
455
456 for (k = 0 ; k < Current.NbrHar ; k += 2) {
457 V->Val[MAX_DIM* k ] = tmp ;
458 V->Val[MAX_DIM*(k+1)] = 0. ;
459 }
460
461 }
462 V->Type = SCALAR ;
463 }
464
465 /* ------------------------------------------------------------------------ */
466 /* Create a Complex Value from k Real Values (of same type!) */
467 /* ------------------------------------------------------------------------ */
468
F_Complex(F_ARG)469 void F_Complex(F_ARG)
470 {
471 int NbrPar = Fct->NbrParameters ;
472 int NbrArg = Fct->NbrArguments ;
473
474 if(NbrArg){
475 if(NbrArg > NBR_MAX_HARMONIC){
476 Message::Error("Too many arguments for Complex[expression-list]{}");
477 return;
478 }
479 }
480 else if(NbrPar){
481 if(NbrPar > NBR_MAX_HARMONIC){
482 Message::Error("Too many parameters for Complex[]{expression-cst-list}");
483 return;
484 }
485 }
486 else{
487 Message::Error("Missing arguments or parameters for Complex[expression-list]{expression-cst-list}");
488 return;
489 }
490
491 int k ;
492
493 if(NbrPar){
494 for (k = 0 ; k < Current.NbrHar ; k++)
495 V->Val[MAX_DIM*k] = Fct->Para[k] ;
496 for (k = Current.NbrHar ; k < NbrPar ; k++)
497 V->Val[MAX_DIM*k] = 0. ;
498 V->Type = SCALAR;
499 return;
500 }
501
502 switch(A->Type){
503
504 case SCALAR :
505 for (k = 0 ; k < Current.NbrHar ; k++) {
506 if((A+k)->Type != A->Type)
507 Message::Error("Mixed type of arguments in function 'Complex'");
508 V->Val[MAX_DIM*k] = (A+k)->Val[0] ;
509 }
510 for (k = Current.NbrHar ; k < NbrArg ; k++) {
511 V->Val[MAX_DIM*k] = 0. ;
512 }
513 break;
514
515 case VECTOR :
516 case TENSOR_DIAG :
517 for (k = 0 ; k < Current.NbrHar ; k++) {
518 if((A+k)->Type != A->Type)
519 Message::Error("Mixed type of arguments in function 'Complex'");
520 V->Val[MAX_DIM*k ] = (A+k)->Val[0] ;
521 V->Val[MAX_DIM*k+1] = (A+k)->Val[1] ;
522 V->Val[MAX_DIM*k+2] = (A+k)->Val[2] ;
523 }
524 for (k = Current.NbrHar ; k < NbrArg ; k++) {
525 V->Val[MAX_DIM*k ] = 0. ;
526 V->Val[MAX_DIM*k+1] = 0. ;
527 V->Val[MAX_DIM*k+2] = 0. ;
528 }
529 break;
530
531 case TENSOR_SYM :
532 for (k = 0 ; k < Current.NbrHar ; k++) {
533 if((A+k)->Type != A->Type)
534 Message::Error("Mixed type of arguments in function 'Complex'");
535 V->Val[MAX_DIM*k ] = (A+k)->Val[0] ;
536 V->Val[MAX_DIM*k+1] = (A+k)->Val[1] ;
537 V->Val[MAX_DIM*k+2] = (A+k)->Val[2] ;
538 V->Val[MAX_DIM*k+3] = (A+k)->Val[3] ;
539 V->Val[MAX_DIM*k+4] = (A+k)->Val[4] ;
540 V->Val[MAX_DIM*k+5] = (A+k)->Val[5] ;
541 }
542 for (k = Current.NbrHar ; k < NbrArg ; k++) {
543 V->Val[MAX_DIM*k ] = 0. ;
544 V->Val[MAX_DIM*k+1] = 0. ;
545 V->Val[MAX_DIM*k+2] = 0. ;
546 V->Val[MAX_DIM*k+3] = 0. ;
547 V->Val[MAX_DIM*k+4] = 0. ;
548 V->Val[MAX_DIM*k+5] = 0. ;
549 }
550 break;
551
552 case TENSOR :
553 for (k = 0 ; k < Current.NbrHar ; k++) {
554 if((A+k)->Type != A->Type)
555 Message::Error("Mixed type of arguments in function 'Complex'");
556 V->Val[MAX_DIM*k ] = (A+k)->Val[0] ;
557 V->Val[MAX_DIM*k+1] = (A+k)->Val[1] ;
558 V->Val[MAX_DIM*k+2] = (A+k)->Val[2] ;
559 V->Val[MAX_DIM*k+3] = (A+k)->Val[3] ;
560 V->Val[MAX_DIM*k+4] = (A+k)->Val[4] ;
561 V->Val[MAX_DIM*k+5] = (A+k)->Val[5] ;
562 V->Val[MAX_DIM*k+6] = (A+k)->Val[6] ;
563 V->Val[MAX_DIM*k+7] = (A+k)->Val[7] ;
564 V->Val[MAX_DIM*k+8] = (A+k)->Val[8] ;
565 }
566 for (k = Current.NbrHar ; k < NbrArg ; k++) {
567 V->Val[MAX_DIM*k ] = 0. ;
568 V->Val[MAX_DIM*k+1] = 0. ;
569 V->Val[MAX_DIM*k+2] = 0. ;
570 V->Val[MAX_DIM*k+3] = 0. ;
571 V->Val[MAX_DIM*k+4] = 0. ;
572 V->Val[MAX_DIM*k+5] = 0. ;
573 V->Val[MAX_DIM*k+6] = 0. ;
574 V->Val[MAX_DIM*k+7] = 0. ;
575 V->Val[MAX_DIM*k+8] = 0. ;
576 }
577 break;
578
579 default :
580 Message::Error("Unknown type of arguments in function 'Complex'");
581 break;
582 }
583
584 V->Type = A->Type ;
585 }
586
587
588 /* ----------------------------------------------------------------------- */
589 /* Get the Real Part of a Value */
590 /* ------------------------------------------------------------------------ */
591
F_Re(F_ARG)592 void F_Re(F_ARG)
593 {
594 int k ;
595
596 switch (A->Type) {
597 case SCALAR :
598 for (k = 0 ; k < Current.NbrHar ; k+=2) {
599 V->Val[MAX_DIM*k] = A->Val[MAX_DIM*k] ;
600 V->Val[MAX_DIM*(k+1)] = 0. ;
601 }
602 break;
603
604 case VECTOR :
605 case TENSOR_DIAG :
606 for (k = 0 ; k < Current.NbrHar ; k+=2) {
607 V->Val[MAX_DIM*k ] = A->Val[MAX_DIM*k ] ;
608 V->Val[MAX_DIM*k+1] = A->Val[MAX_DIM*k+1] ;
609 V->Val[MAX_DIM*k+2] = A->Val[MAX_DIM*k+2] ;
610 V->Val[MAX_DIM*(k+1) ] = 0. ;
611 V->Val[MAX_DIM*(k+1)+1] = 0. ;
612 V->Val[MAX_DIM*(k+1)+2] = 0. ;
613 }
614 break;
615
616 case TENSOR_SYM :
617 for (k = 0 ; k < Current.NbrHar ; k+=2) {
618 V->Val[MAX_DIM*k ] = A->Val[MAX_DIM*k ] ;
619 V->Val[MAX_DIM*k+1] = A->Val[MAX_DIM*k+1] ;
620 V->Val[MAX_DIM*k+2] = A->Val[MAX_DIM*k+2] ;
621 V->Val[MAX_DIM*k+3] = A->Val[MAX_DIM*k+3] ;
622 V->Val[MAX_DIM*k+4] = A->Val[MAX_DIM*k+4] ;
623 V->Val[MAX_DIM*k+5] = A->Val[MAX_DIM*k+5] ;
624 V->Val[MAX_DIM*(k+1) ] = 0. ;
625 V->Val[MAX_DIM*(k+1)+1] = 0. ;
626 V->Val[MAX_DIM*(k+1)+2] = 0. ;
627 V->Val[MAX_DIM*(k+1)+3] = 0. ;
628 V->Val[MAX_DIM*(k+1)+4] = 0. ;
629 V->Val[MAX_DIM*(k+1)+5] = 0. ;
630 }
631 break;
632
633 case TENSOR :
634 for (k = 0 ; k < Current.NbrHar ; k+=2) {
635 V->Val[MAX_DIM*k ] = A->Val[MAX_DIM*k ] ;
636 V->Val[MAX_DIM*k+1] = A->Val[MAX_DIM*k+1] ;
637 V->Val[MAX_DIM*k+2] = A->Val[MAX_DIM*k+2] ;
638 V->Val[MAX_DIM*k+3] = A->Val[MAX_DIM*k+3] ;
639 V->Val[MAX_DIM*k+4] = A->Val[MAX_DIM*k+4] ;
640 V->Val[MAX_DIM*k+5] = A->Val[MAX_DIM*k+5] ;
641 V->Val[MAX_DIM*k+6] = A->Val[MAX_DIM*k+6] ;
642 V->Val[MAX_DIM*k+7] = A->Val[MAX_DIM*k+7] ;
643 V->Val[MAX_DIM*k+8] = A->Val[MAX_DIM*k+8] ;
644 V->Val[MAX_DIM*(k+1) ] = 0. ;
645 V->Val[MAX_DIM*(k+1)+1] = 0. ;
646 V->Val[MAX_DIM*(k+1)+2] = 0. ;
647 V->Val[MAX_DIM*(k+1)+3] = 0. ;
648 V->Val[MAX_DIM*(k+1)+4] = 0. ;
649 V->Val[MAX_DIM*(k+1)+5] = 0. ;
650 V->Val[MAX_DIM*(k+1)+6] = 0. ;
651 V->Val[MAX_DIM*(k+1)+7] = 0. ;
652 V->Val[MAX_DIM*(k+1)+8] = 0. ;
653 }
654 break;
655
656 default :
657 Message::Error("Unknown type of arguments in function 'Re'");
658 break;
659 }
660
661 V->Type = A->Type ;
662 }
663
664 /* ------------------------------------------------------------------------ */
665 /* Get the Imaginary Part of a Value */
666 /* ------------------------------------------------------------------------ */
667
F_Im(F_ARG)668 void F_Im(F_ARG)
669 {
670 int k ;
671
672 switch (A->Type) {
673 case SCALAR :
674 for (k = 0 ; k < Current.NbrHar ; k+=2) {
675 V->Val[MAX_DIM*k] = A->Val[MAX_DIM*(k+1)] ;
676 V->Val[MAX_DIM*(k+1)] = 0. ;
677 }
678 break;
679
680 case VECTOR :
681 case TENSOR_DIAG :
682 for (k = 0 ; k < Current.NbrHar ; k+=2) {
683 V->Val[MAX_DIM*k ] = A->Val[MAX_DIM*(k+1) ] ;
684 V->Val[MAX_DIM*k+1] = A->Val[MAX_DIM*(k+1)+1] ;
685 V->Val[MAX_DIM*k+2] = A->Val[MAX_DIM*(k+1)+2] ;
686 V->Val[MAX_DIM*(k+1) ] = 0. ;
687 V->Val[MAX_DIM*(k+1)+1] = 0. ;
688 V->Val[MAX_DIM*(k+1)+2] = 0. ;
689 }
690 break;
691
692 case TENSOR_SYM :
693 for (k = 0 ; k < Current.NbrHar ; k+=2) {
694 V->Val[MAX_DIM*k ] = A->Val[MAX_DIM*(k+1) ] ;
695 V->Val[MAX_DIM*k+1] = A->Val[MAX_DIM*(k+1)+1] ;
696 V->Val[MAX_DIM*k+2] = A->Val[MAX_DIM*(k+1)+2] ;
697 V->Val[MAX_DIM*k+3] = A->Val[MAX_DIM*(k+1)+3] ;
698 V->Val[MAX_DIM*k+4] = A->Val[MAX_DIM*(k+1)+4] ;
699 V->Val[MAX_DIM*k+5] = A->Val[MAX_DIM*(k+1)+5] ;
700 V->Val[MAX_DIM*(k+1) ] = 0. ;
701 V->Val[MAX_DIM*(k+1)+1] = 0. ;
702 V->Val[MAX_DIM*(k+1)+2] = 0. ;
703 V->Val[MAX_DIM*(k+1)+3] = 0. ;
704 V->Val[MAX_DIM*(k+1)+4] = 0. ;
705 V->Val[MAX_DIM*(k+1)+5] = 0. ;
706 }
707 break;
708
709 case TENSOR :
710 for (k = 0 ; k < Current.NbrHar ; k+=2) {
711 V->Val[MAX_DIM*k ] = A->Val[MAX_DIM*(k+1) ] ;
712 V->Val[MAX_DIM*k+1] = A->Val[MAX_DIM*(k+1)+1] ;
713 V->Val[MAX_DIM*k+2] = A->Val[MAX_DIM*(k+1)+2] ;
714 V->Val[MAX_DIM*k+3] = A->Val[MAX_DIM*(k+1)+3] ;
715 V->Val[MAX_DIM*k+4] = A->Val[MAX_DIM*(k+1)+4] ;
716 V->Val[MAX_DIM*k+5] = A->Val[MAX_DIM*(k+1)+5] ;
717 V->Val[MAX_DIM*k+6] = A->Val[MAX_DIM*(k+1)+6] ;
718 V->Val[MAX_DIM*k+7] = A->Val[MAX_DIM*(k+1)+7] ;
719 V->Val[MAX_DIM*k+8] = A->Val[MAX_DIM*(k+1)+8] ;
720 V->Val[MAX_DIM*(k+1) ] = 0. ;
721 V->Val[MAX_DIM*(k+1)+1] = 0. ;
722 V->Val[MAX_DIM*(k+1)+2] = 0. ;
723 V->Val[MAX_DIM*(k+1)+3] = 0. ;
724 V->Val[MAX_DIM*(k+1)+4] = 0. ;
725 V->Val[MAX_DIM*(k+1)+5] = 0. ;
726 V->Val[MAX_DIM*(k+1)+6] = 0. ;
727 V->Val[MAX_DIM*(k+1)+7] = 0. ;
728 V->Val[MAX_DIM*(k+1)+8] = 0. ;
729 }
730 break;
731
732 default :
733 Message::Error("Unknown type of arguments in function 'Im'");
734 break;
735 }
736
737 V->Type = A->Type ;
738 }
739
740 /* ------------------------------------------------------------------------ */
741 /* Conjugate */
742 /* ------------------------------------------------------------------------ */
743
F_Conj(F_ARG)744 void F_Conj(F_ARG)
745 {
746 int k ;
747
748 switch (A->Type) {
749 case SCALAR :
750 for (k = 0 ; k < Current.NbrHar ; k+=2) {
751 V->Val[MAX_DIM*k] = A->Val[MAX_DIM*k] ;
752 V->Val[MAX_DIM*(k+1)] = -A->Val[MAX_DIM*(k+1)] ;
753 }
754 break;
755
756 case VECTOR :
757 case TENSOR_DIAG :
758 for (k = 0 ; k < Current.NbrHar ; k+=2) {
759 V->Val[MAX_DIM*k ] = A->Val[MAX_DIM*k ] ;
760 V->Val[MAX_DIM*k+1] = A->Val[MAX_DIM*k+1] ;
761 V->Val[MAX_DIM*k+2] = A->Val[MAX_DIM*k+2] ;
762 V->Val[MAX_DIM*(k+1) ] = -A->Val[MAX_DIM*(k+1) ] ;
763 V->Val[MAX_DIM*(k+1)+1] = -A->Val[MAX_DIM*(k+1)+1] ;
764 V->Val[MAX_DIM*(k+1)+2] = -A->Val[MAX_DIM*(k+1)+2] ;
765 }
766 break;
767
768 case TENSOR_SYM :
769 for (k = 0 ; k < Current.NbrHar ; k+=2) {
770 V->Val[MAX_DIM*k ] = A->Val[MAX_DIM*k ] ;
771 V->Val[MAX_DIM*k+1] = A->Val[MAX_DIM*k+1] ;
772 V->Val[MAX_DIM*k+2] = A->Val[MAX_DIM*k+2] ;
773 V->Val[MAX_DIM*k+3] = A->Val[MAX_DIM*k+3] ;
774 V->Val[MAX_DIM*k+4] = A->Val[MAX_DIM*k+4] ;
775 V->Val[MAX_DIM*k+5] = A->Val[MAX_DIM*k+5] ;
776 V->Val[MAX_DIM*(k+1) ] = -A->Val[MAX_DIM*(k+1) ] ;
777 V->Val[MAX_DIM*(k+1)+1] = -A->Val[MAX_DIM*(k+1)+1] ;
778 V->Val[MAX_DIM*(k+1)+2] = -A->Val[MAX_DIM*(k+1)+2] ;
779 V->Val[MAX_DIM*(k+1)+3] = -A->Val[MAX_DIM*(k+1)+3] ;
780 V->Val[MAX_DIM*(k+1)+4] = -A->Val[MAX_DIM*(k+1)+4] ;
781 V->Val[MAX_DIM*(k+1)+5] = -A->Val[MAX_DIM*(k+1)+5] ;
782 }
783 break;
784
785 case TENSOR :
786 for (k = 0 ; k < Current.NbrHar ; k+=2) {
787 V->Val[MAX_DIM*k ] = A->Val[MAX_DIM*k ] ;
788 V->Val[MAX_DIM*k+1] = A->Val[MAX_DIM*k+1] ;
789 V->Val[MAX_DIM*k+2] = A->Val[MAX_DIM*k+2] ;
790 V->Val[MAX_DIM*k+3] = A->Val[MAX_DIM*k+3] ;
791 V->Val[MAX_DIM*k+4] = A->Val[MAX_DIM*k+4] ;
792 V->Val[MAX_DIM*k+5] = A->Val[MAX_DIM*k+5] ;
793 V->Val[MAX_DIM*k+6] = A->Val[MAX_DIM*k+6] ;
794 V->Val[MAX_DIM*k+7] = A->Val[MAX_DIM*k+7] ;
795 V->Val[MAX_DIM*k+8] = A->Val[MAX_DIM*k+8] ;
796 V->Val[MAX_DIM*(k+1) ] = -A->Val[MAX_DIM*(k+1) ] ;
797 V->Val[MAX_DIM*(k+1)+1] = -A->Val[MAX_DIM*(k+1)+1] ;
798 V->Val[MAX_DIM*(k+1)+2] = -A->Val[MAX_DIM*(k+1)+2] ;
799 V->Val[MAX_DIM*(k+1)+3] = -A->Val[MAX_DIM*(k+1)+3] ;
800 V->Val[MAX_DIM*(k+1)+4] = -A->Val[MAX_DIM*(k+1)+4] ;
801 V->Val[MAX_DIM*(k+1)+5] = -A->Val[MAX_DIM*(k+1)+5] ;
802 V->Val[MAX_DIM*(k+1)+6] = -A->Val[MAX_DIM*(k+1)+6] ;
803 V->Val[MAX_DIM*(k+1)+7] = -A->Val[MAX_DIM*(k+1)+7] ;
804 V->Val[MAX_DIM*(k+1)+8] = -A->Val[MAX_DIM*(k+1)+8] ;
805 }
806 break;
807
808 default :
809 Message::Error("Unknown type of arguments in function 'Conj'");
810 break;
811 }
812
813 V->Type = A->Type ;
814 }
815
816 /* -------------------------------------------------------------------------------- */
817 /* Cartesian coordinates (Re,Im) to polar coordinates (Amplitude,phase[Radians]) */
818 /* -------------------------------------------------------------------------------- */
819
F_Cart2Pol(F_ARG)820 void F_Cart2Pol(F_ARG)
821 {
822 int k ;
823 double Re, Im;
824
825 switch (A->Type) {
826 case SCALAR :
827 for (k = 0 ; k < Current.NbrHar ; k+=2) {
828 Re = A->Val[MAX_DIM*k] ; Im = A->Val[MAX_DIM*(k+1)] ;
829 V->Val[MAX_DIM*k] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)] = atan2(Im,Re);
830 }
831 break;
832
833 case VECTOR :
834 case TENSOR_DIAG :
835 for (k = 0 ; k < Current.NbrHar ; k+=2) {
836 Re = A->Val[MAX_DIM*k ] ; Im = A->Val[MAX_DIM*(k+1) ] ;
837 V->Val[MAX_DIM*k ] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1) ] = atan2(Im,Re);
838 Re = A->Val[MAX_DIM*k+1] ; Im = A->Val[MAX_DIM*(k+1)+1] ;
839 V->Val[MAX_DIM*k+1] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+1] = atan2(Im,Re);
840 Re = A->Val[MAX_DIM*k+2] ; Im = A->Val[MAX_DIM*(k+1)+2] ;
841 V->Val[MAX_DIM*k+2] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+2] = atan2(Im,Re);
842 }
843 break;
844
845 case TENSOR_SYM :
846 for (k = 0 ; k < Current.NbrHar ; k+=2) {
847 Re = A->Val[MAX_DIM*k ] ; Im = A->Val[MAX_DIM*(k+1) ] ;
848 V->Val[MAX_DIM*k ] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1) ] = atan2(Im,Re);
849 Re = A->Val[MAX_DIM*k+1] ; Im = A->Val[MAX_DIM*(k+1)+1] ;
850 V->Val[MAX_DIM*k+1] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+1] = atan2(Im,Re);
851 Re = A->Val[MAX_DIM*k+2] ; Im = A->Val[MAX_DIM*(k+1)+2] ;
852 V->Val[MAX_DIM*k+2] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+2] = atan2(Im,Re);
853 Re = A->Val[MAX_DIM*k+3] ; Im = A->Val[MAX_DIM*(k+1)+3] ;
854 V->Val[MAX_DIM*k+3] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+3] = atan2(Im,Re);
855 Re = A->Val[MAX_DIM*k+4] ; Im = A->Val[MAX_DIM*(k+1)+4] ;
856 V->Val[MAX_DIM*k+4] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+4] = atan2(Im,Re);
857 Re = A->Val[MAX_DIM*k+5] ; Im = A->Val[MAX_DIM*(k+1)+5] ;
858 V->Val[MAX_DIM*k+5] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+5] = atan2(Im,Re);
859 }
860 break;
861
862 case TENSOR :
863 for (k = 0 ; k < Current.NbrHar ; k+=2) {
864 Re = A->Val[MAX_DIM*k ] ; Im = A->Val[MAX_DIM*(k+1) ] ;
865 V->Val[MAX_DIM*k ] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1) ] = atan2(Im,Re);
866 Re = A->Val[MAX_DIM*k+1] ; Im = A->Val[MAX_DIM*(k+1)+1] ;
867 V->Val[MAX_DIM*k+1] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+1] = atan2(Im,Re);
868 Re = A->Val[MAX_DIM*k+2] ; Im = A->Val[MAX_DIM*(k+1)+2] ;
869 V->Val[MAX_DIM*k+2] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+2] = atan2(Im,Re);
870 Re = A->Val[MAX_DIM*k+3] ; Im = A->Val[MAX_DIM*(k+1)+3] ;
871 V->Val[MAX_DIM*k+3] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+3] = atan2(Im,Re);
872 Re = A->Val[MAX_DIM*k+4] ; Im = A->Val[MAX_DIM*(k+1)+4] ;
873 V->Val[MAX_DIM*k+4] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+4] = atan2(Im,Re);
874 Re = A->Val[MAX_DIM*k+5] ; Im = A->Val[MAX_DIM*(k+1)+5] ;
875 V->Val[MAX_DIM*k+5] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+5] = atan2(Im,Re);
876 Re = A->Val[MAX_DIM*k+6] ; Im = A->Val[MAX_DIM*(k+1)+6] ;
877 V->Val[MAX_DIM*k+6] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+6] = atan2(Im,Re);
878 Re = A->Val[MAX_DIM*k+7] ; Im = A->Val[MAX_DIM*(k+1)+7] ;
879 V->Val[MAX_DIM*k+7] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+7] = atan2(Im,Re);
880 Re = A->Val[MAX_DIM*k+8] ; Im = A->Val[MAX_DIM*(k+1)+8] ;
881 V->Val[MAX_DIM*k+8] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+8] = atan2(Im,Re);
882 }
883 break;
884
885 default :
886 Message::Error("Unknown type of arguments in function 'Cart2Pol'");
887 break;
888 }
889
890 V->Type = A->Type ;
891 }
892
893 /* ------------------------------------------------------------------------ */
894 /* Create 1 Vector from 3 Scalar */
895 /* ------------------------------------------------------------------------ */
896
F_Vector(F_ARG)897 void F_Vector(F_ARG)
898 {
899 int k ;
900
901 if(A->Type != SCALAR || (A+1)->Type != SCALAR || (A+2)->Type != SCALAR)
902 Message::Error("Non scalar argument(s) for function 'Vector'");
903
904 for (k = 0 ; k < Current.NbrHar ; k++) {
905 V->Val[MAX_DIM*k ] = (A )->Val[MAX_DIM*k] ;
906 V->Val[MAX_DIM*k+1] = (A+1)->Val[MAX_DIM*k] ;
907 V->Val[MAX_DIM*k+2] = (A+2)->Val[MAX_DIM*k] ;
908 }
909 V->Type = VECTOR ;
910 }
911
912 /* ------------------------------------------------------------------------ */
913 /* Create 1 Tensor from 9 Scalar */
914 /* ------------------------------------------------------------------------ */
915
F_Tensor(F_ARG)916 void F_Tensor(F_ARG)
917 {
918 int k ;
919
920 if( (A)->Type != SCALAR || (A+1)->Type != SCALAR || (A+2)->Type != SCALAR ||
921 (A+3)->Type != SCALAR || (A+4)->Type != SCALAR || (A+5)->Type != SCALAR ||
922 (A+6)->Type != SCALAR || (A+7)->Type != SCALAR || (A+8)->Type != SCALAR )
923 Message::Error("Non scalar argument(s) for function 'Tensor'");
924
925 for (k = 0 ; k < Current.NbrHar ; k++) {
926 V->Val[MAX_DIM*k ] = (A )->Val[MAX_DIM*k] ;
927 V->Val[MAX_DIM*k+1] = (A+1)->Val[MAX_DIM*k] ;
928 V->Val[MAX_DIM*k+2] = (A+2)->Val[MAX_DIM*k] ;
929 V->Val[MAX_DIM*k+3] = (A+3)->Val[MAX_DIM*k] ;
930 V->Val[MAX_DIM*k+4] = (A+4)->Val[MAX_DIM*k] ;
931 V->Val[MAX_DIM*k+5] = (A+5)->Val[MAX_DIM*k] ;
932 V->Val[MAX_DIM*k+6] = (A+6)->Val[MAX_DIM*k] ;
933 V->Val[MAX_DIM*k+7] = (A+7)->Val[MAX_DIM*k] ;
934 V->Val[MAX_DIM*k+8] = (A+8)->Val[MAX_DIM*k] ;
935 }
936 V->Type = TENSOR ;
937 }
938
939 /* ------------------------------------------------------------------------ */
940 /* Create 1 Symmetric Tensor from 6 Scalar */
941 /* ------------------------------------------------------------------------ */
942
F_TensorSym(F_ARG)943 void F_TensorSym(F_ARG)
944 {
945 int k ;
946
947 if( (A)->Type != SCALAR || (A+1)->Type != SCALAR || (A+2)->Type != SCALAR ||
948 (A+3)->Type != SCALAR || (A+4)->Type != SCALAR || (A+5)->Type != SCALAR )
949 Message::Error("Non scalar argument(s) for function 'TensorSym'");
950
951 for (k = 0 ; k < Current.NbrHar ; k++) {
952 V->Val[MAX_DIM*k ] = (A )->Val[MAX_DIM*k] ;
953 V->Val[MAX_DIM*k+1] = (A+1)->Val[MAX_DIM*k] ;
954 V->Val[MAX_DIM*k+2] = (A+2)->Val[MAX_DIM*k] ;
955 V->Val[MAX_DIM*k+3] = (A+3)->Val[MAX_DIM*k] ;
956 V->Val[MAX_DIM*k+4] = (A+4)->Val[MAX_DIM*k] ;
957 V->Val[MAX_DIM*k+5] = (A+5)->Val[MAX_DIM*k] ;
958 }
959 V->Type = TENSOR_SYM ;
960 }
961
962 /* ------------------------------------------------------------------------ */
963 /* Create 1 Diagonal Tensor from 3 Scalar */
964 /* ------------------------------------------------------------------------ */
965
F_TensorDiag(F_ARG)966 void F_TensorDiag(F_ARG)
967 {
968 int k ;
969
970 if(A->Type != SCALAR || (A+1)->Type != SCALAR || (A+2)->Type != SCALAR)
971 Message::Error("Non scalar argument(s) for function 'TensorDiag'");
972
973 for (k = 0 ; k < Current.NbrHar ; k++) {
974 V->Val[MAX_DIM*k ] = A->Val[MAX_DIM*k] ;
975 V->Val[MAX_DIM*k+1] = (A+1)->Val[MAX_DIM*k] ;
976 V->Val[MAX_DIM*k+2] = (A+2)->Val[MAX_DIM*k] ;
977 }
978 V->Type = TENSOR_DIAG ;
979 }
980
981 /* ------------------------------------------------------------------------ */
982 /* Create 1 Tensor from 3 Vector */
983 /* ------------------------------------------------------------------------ */
984
F_TensorV(F_ARG)985 void F_TensorV(F_ARG)
986 {
987 int k ;
988
989 if((A)->Type != VECTOR || (A+1)->Type != VECTOR || (A+2)->Type != VECTOR)
990 Message::Error("Non scalar argument(s) for function 'TensorV'");
991
992 for (k = 0 ; k < Current.NbrHar ; k++) {
993 V->Val[MAX_DIM*k ] = (A )->Val[MAX_DIM*k ] ;
994 V->Val[MAX_DIM*k+1] = (A )->Val[MAX_DIM*k+1] ;
995 V->Val[MAX_DIM*k+2] = (A )->Val[MAX_DIM*k+2] ;
996 V->Val[MAX_DIM*k+3] = (A+1)->Val[MAX_DIM*k ] ;
997 V->Val[MAX_DIM*k+4] = (A+1)->Val[MAX_DIM*k+1] ;
998 V->Val[MAX_DIM*k+5] = (A+1)->Val[MAX_DIM*k+2] ;
999 V->Val[MAX_DIM*k+6] = (A+2)->Val[MAX_DIM*k ] ;
1000 V->Val[MAX_DIM*k+7] = (A+2)->Val[MAX_DIM*k+1] ;
1001 V->Val[MAX_DIM*k+8] = (A+2)->Val[MAX_DIM*k+2] ;
1002 }
1003 V->Type = TENSOR ;
1004 }
1005
1006 /* ------------------------------------------------------------------------ */
1007 /* Dyadic product */
1008 /* ------------------------------------------------------------------------ */
1009
F_SquDyadicProduct(F_ARG)1010 void F_SquDyadicProduct(F_ARG)
1011 {
1012 int k ;
1013 double t11, t12, t13, t22, t23, t33 ;
1014
1015 if (A->Type != VECTOR)
1016 Message::Error("Non vector argument for function 'TensorDyadic'");
1017
1018 t11 = SQU(A->Val[0]) ; t22 = SQU(A->Val[1]) ; t33 = SQU(A->Val[2]) ;
1019 t12 = A->Val[0] * A->Val[1] ; t13 = A->Val[0] * A->Val[2] ;
1020 t23 = A->Val[1] * A->Val[2] ;
1021
1022 V->Val[0] = t11 ; V->Val[1] = t12 ; V->Val[2] = t13 ;
1023 V->Val[3] = t22 ; V->Val[4] = t23 ; V->Val[5] = t33 ;
1024
1025 /* Attention : a revoir */
1026 if (Current.NbrHar > 1) {
1027 V->Val[MAX_DIM ] = V->Val[MAX_DIM+1] = V->Val[MAX_DIM+2] =
1028 V->Val[MAX_DIM+3] = V->Val[MAX_DIM+4] = V->Val[MAX_DIM+5] = 0. ;
1029 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k++) {
1030 V->Val[MAX_DIM*k ] = V->Val[MAX_DIM*k+1] = V->Val[MAX_DIM*k+2] =
1031 V->Val[MAX_DIM*k+3] = V->Val[MAX_DIM*k+4] = V->Val[MAX_DIM*k+5] = 0. ;
1032 }
1033 }
1034
1035 V->Type = TENSOR_SYM ;
1036 }
1037
1038 /* ------------------------------------------------------------------------ */
1039 /* Get Vector Components */
1040 /* ------------------------------------------------------------------------ */
1041
1042 #define get_comp_vector(index, string) \
1043 int k ; \
1044 \
1045 if(A->Type != VECTOR) \
1046 Message::Error("Non vector argument for function '" string "'"); \
1047 \
1048 for (k = 0 ; k < Current.NbrHar ; k++) { \
1049 V->Val[MAX_DIM*k ] = A->Val[MAX_DIM*k+index] ; \
1050 } \
1051 V->Type = SCALAR ;
1052
F_CompX(F_ARG)1053 void F_CompX(F_ARG){ get_comp_vector(0, "CompX") }
F_CompY(F_ARG)1054 void F_CompY(F_ARG){ get_comp_vector(1, "CompY") }
F_CompZ(F_ARG)1055 void F_CompZ(F_ARG){ get_comp_vector(2, "CompZ") }
1056
F_Comp(F_ARG)1057 void F_Comp(F_ARG){
1058 if (Fct->NbrParameters != 1)
1059 Message::Error("Function 'Comp': one parameter needed to define component index");
1060 if ((int)(Fct->Para[0]) < 0 || (int)(Fct->Para[0]) > 2)
1061 Message::Error("Function 'Comp': parameter (%g) out of range (must be 0, 1 or 2)", Fct->Para[0]);
1062
1063 get_comp_vector((int)(Fct->Para[0]), "Comp")
1064 }
1065
1066 #undef get_comp_vector
1067
1068 /* ------------------------------------------------------------------------ */
1069 /* Get Tensor Components */
1070 /* ------------------------------------------------------------------------ */
1071
1072 #define get_comp_tensor(i, is, id, string) \
1073 int k ; \
1074 \
1075 switch(A->Type) { \
1076 case TENSOR : \
1077 for (k=0; k<Current.NbrHar; k++) V->Val[MAX_DIM*k] = A->Val[MAX_DIM*k+(i)] ; \
1078 break ; \
1079 case TENSOR_SYM : \
1080 for (k=0; k<Current.NbrHar; k++) V->Val[MAX_DIM*k] = A->Val[MAX_DIM*k+(is)] ; \
1081 break ; \
1082 case TENSOR_DIAG : \
1083 if(id >= 0) \
1084 for (k=0; k<Current.NbrHar; k++) V->Val[MAX_DIM*k] = A->Val[MAX_DIM*k+(id)] ; \
1085 else \
1086 for (k=0; k<Current.NbrHar; k++) V->Val[MAX_DIM*k] = 0.; \
1087 break ; \
1088 case SCALAR : \
1089 for (k=0; k<Current.NbrHar; k++) V->Val[MAX_DIM*k] = A->Val[MAX_DIM*k] ; \
1090 break ; \
1091 default : \
1092 Message::Error("Non tensor or scalar argument for function '" string "'"); \
1093 break; \
1094 } \
1095 V->Type = SCALAR ;
1096
F_CompXX(F_ARG)1097 void F_CompXX(F_ARG){ get_comp_tensor(0,0, 0,"CompXX") }
F_CompXY(F_ARG)1098 void F_CompXY(F_ARG){ get_comp_tensor(1,1,-1,"CompXY") }
F_CompXZ(F_ARG)1099 void F_CompXZ(F_ARG){ get_comp_tensor(2,2,-1,"CompXZ") }
F_CompYX(F_ARG)1100 void F_CompYX(F_ARG){ get_comp_tensor(3,1,-1,"CompYX") }
F_CompYY(F_ARG)1101 void F_CompYY(F_ARG){ get_comp_tensor(4,3, 1,"CompYY") }
F_CompYZ(F_ARG)1102 void F_CompYZ(F_ARG){ get_comp_tensor(5,4,-1,"CompYZ") }
F_CompZX(F_ARG)1103 void F_CompZX(F_ARG){ get_comp_tensor(6,2,-1,"CompZX") }
F_CompZY(F_ARG)1104 void F_CompZY(F_ARG){ get_comp_tensor(7,4,-1,"CompZY") }
F_CompZZ(F_ARG)1105 void F_CompZZ(F_ARG){ get_comp_tensor(8,5, 2,"CompZZ") }
1106
1107 #undef get_comp_tensor
1108
1109 /* ------------------------------------------------------------------------ */
1110 /* Get Tensor for transformation of vector */
1111 /* from cartesian to spherical coordinate system */
1112 /* ------------------------------------------------------------------------ */
1113
F_Cart2Sph(F_ARG)1114 void F_Cart2Sph(F_ARG)
1115 {
1116 int k ;
1117 double theta, phi ;
1118
1119 if((A)->Type != VECTOR)
1120 Message::Error("Vector argument required for Function 'Cart2Sph'");
1121
1122 /* Warning! This is the physic's convention. For the math
1123 convention, switch theta and phi. */
1124
1125 theta = atan2( sqrt(SQU(A->Val[0])+ SQU(A->Val[1])) , A->Val[2] ) ;
1126 phi = atan2( A->Val[1] , A->Val[0] ) ;
1127
1128 /* r basis vector */
1129 V->Val[0] = sin(theta) * cos(phi) ;
1130 V->Val[1] = sin(theta) * sin(phi) ;
1131 V->Val[2] = cos(theta) ;
1132
1133 /* theta basis vector */
1134 V->Val[3] = cos(theta) * cos(phi) ;
1135 V->Val[4] = cos(theta) * sin(phi) ;
1136 V->Val[5] = - sin(theta) ;
1137
1138 /* phi basis vector */
1139 V->Val[6] = - sin(phi) ;
1140 V->Val[7] = cos(phi) ;
1141 V->Val[8] = 0. ;
1142
1143 for (k = 0 ; k < Current.NbrHar ; k+=2) {
1144 V->Val[MAX_DIM*k ] = V->Val[0] ;
1145 V->Val[MAX_DIM*k+1] = V->Val[1] ;
1146 V->Val[MAX_DIM*k+2] = V->Val[2] ;
1147 V->Val[MAX_DIM*k+3] = V->Val[3] ;
1148 V->Val[MAX_DIM*k+4] = V->Val[4] ;
1149 V->Val[MAX_DIM*k+5] = V->Val[5] ;
1150 V->Val[MAX_DIM*k+6] = V->Val[6] ;
1151 V->Val[MAX_DIM*k+7] = V->Val[7] ;
1152 V->Val[MAX_DIM*k+8] = V->Val[8] ;
1153 V->Val[MAX_DIM*(k+1) ] = 0. ;
1154 V->Val[MAX_DIM*(k+1)+1] = 0. ;
1155 V->Val[MAX_DIM*(k+1)+2] = 0. ;
1156 V->Val[MAX_DIM*(k+1)+3] = 0. ;
1157 V->Val[MAX_DIM*(k+1)+4] = 0. ;
1158 V->Val[MAX_DIM*(k+1)+5] = 0. ;
1159 V->Val[MAX_DIM*(k+1)+6] = 0. ;
1160 V->Val[MAX_DIM*(k+1)+7] = 0. ;
1161 V->Val[MAX_DIM*(k+1)+8] = 0. ;
1162 }
1163 V->Type = TENSOR ;
1164 }
1165
1166 /* ------------------------------------------------------------------------ */
1167 /* Get Tensor for transformation of vector */
1168 /* from cartesian to cylindric coordinate system */
1169 /* vector -> Cart2Cyl[XYZ[]] * vector */
1170 /* (x,y,z)-components -> (radial, tangential, axial)-components */
1171 /* ------------------------------------------------------------------------ */
1172
F_Cart2Cyl(F_ARG)1173 void F_Cart2Cyl(F_ARG)
1174 {
1175 int k ;
1176 double theta ;
1177
1178 if((A)->Type != VECTOR)
1179 Message::Error("Vector argument required for Function 'Cart2Cyl'");
1180
1181 theta = atan2(A->Val[1] , A->Val[0]) ;
1182
1183 V->Val[0] = cos(theta) ;
1184 V->Val[1] = sin(theta) ;
1185 V->Val[2] = 0 ;
1186 V->Val[3] = -sin(theta) ;
1187 V->Val[4] = cos(theta) ;
1188 V->Val[5] = 0 ;
1189 V->Val[6] = 0 ;
1190 V->Val[7] = 0 ;
1191 V->Val[8] = 1. ;
1192
1193 for (k = 0 ; k < Current.NbrHar ; k+=2) {
1194 V->Val[MAX_DIM*k ] = V->Val[0] ;
1195 V->Val[MAX_DIM*k+1] = V->Val[1] ;
1196 V->Val[MAX_DIM*k+2] = V->Val[2] ;
1197 V->Val[MAX_DIM*k+3] = V->Val[3] ;
1198 V->Val[MAX_DIM*k+4] = V->Val[4] ;
1199 V->Val[MAX_DIM*k+5] = V->Val[5] ;
1200 V->Val[MAX_DIM*k+6] = V->Val[6] ;
1201 V->Val[MAX_DIM*k+7] = V->Val[7] ;
1202 V->Val[MAX_DIM*k+8] = V->Val[8] ;
1203 V->Val[MAX_DIM*(k+1) ] = 0. ;
1204 V->Val[MAX_DIM*(k+1)+1] = 0. ;
1205 V->Val[MAX_DIM*(k+1)+2] = 0. ;
1206 V->Val[MAX_DIM*(k+1)+3] = 0. ;
1207 V->Val[MAX_DIM*(k+1)+4] = 0. ;
1208 V->Val[MAX_DIM*(k+1)+5] = 0. ;
1209 V->Val[MAX_DIM*(k+1)+6] = 0. ;
1210 V->Val[MAX_DIM*(k+1)+7] = 0. ;
1211 V->Val[MAX_DIM*(k+1)+8] = 0. ;
1212 }
1213 V->Type = TENSOR ;
1214 }
1215
1216 /* ------------------------------------------------------------------------ */
1217 /* U n i t V e c t o r X, Y, Z */
1218 /* ------------------------------------------------------------------------ */
1219
F_UnitVectorX(F_ARG)1220 void F_UnitVectorX(F_ARG)
1221 {
1222 int k ;
1223
1224 for (k = 0 ; k < Current.NbrHar ; k++) {
1225 V->Val[MAX_DIM*k ] = (k)? 0.:1. ;
1226 V->Val[MAX_DIM*k+1] = 0. ;
1227 V->Val[MAX_DIM*k+2] = 0. ;
1228 }
1229 V->Type = VECTOR ;
1230 }
1231
F_UnitVectorY(F_ARG)1232 void F_UnitVectorY(F_ARG)
1233 {
1234 int k ;
1235
1236 for (k = 0 ; k < Current.NbrHar ; k++) {
1237 V->Val[MAX_DIM*k ] = 0. ;
1238 V->Val[MAX_DIM*k+1] = (k)? 0.:1. ;
1239 V->Val[MAX_DIM*k+2] = 0. ;
1240 }
1241 V->Type = VECTOR ;
1242 }
1243
F_UnitVectorZ(F_ARG)1244 void F_UnitVectorZ(F_ARG)
1245 {
1246 int k ;
1247
1248 for (k = 0 ; k < Current.NbrHar ; k++) {
1249 V->Val[MAX_DIM*k ] = 0. ;
1250 V->Val[MAX_DIM*k+1] = 0. ;
1251 V->Val[MAX_DIM*k+2] = (k)? 0.:1. ;
1252 }
1253 V->Type = VECTOR ;
1254 }
1255