Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

modnaff.cc @ 3

Last change on this file since 3 was 3, checked in by zhangj, 12 years ago
Initiale import
Property svn:executable set to ``*
File size: 124.6 KB

Line
1	/***************************************************************************
2	modnaff.c - description
3	-------------------
4	begin : Fri Dec 13 20:24:44 CET 2002
5	copyright : (C) 2002 by Laurent Nadolski
6	email : nadolski@synchrotron-soleil.fr
7
8	Origine : Bureau des Longitudes (IMCCE)
9	J. Laskar, M. Gastineau
10
11	***************************************************************************/
12
13
14	/* MODNAFF.C version 0.96 */
15	/* Astronomie et systemes dynamiques */
16	/* J. LASKAR Bureau des Longitudes, Paris (version originale fortran) */
17	/* M. GASTINEAU, portage en C, 03/09/98 */
18	/* v0.96 M. GASTINEAU 09/09/98 : modification dans la fonction naf_ztder */
19	/* dans les boucles de I en I-1. */
20	/* v0.96 M. GASTINEAU 09/11/98 : remplacement des variables PI2 par PICARRE*/
21	/* car PI2 est un define sous LINUX. */
22	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes */
23	/* inlines. */
24	/* v0.96 M. GASTINEAU 07/12/98 : modification dans naf_frefin car bug lors*/
25	/* de l'optimisation sur Mac. */
26	/* v0.96 M. GASTINEAU 18/12/98 : modification de naf_iniwin */
27	/* v0.96 M. GASTINEAU 06/01/99 : ajout de la liberation de la liste des */
28	/* fenetres dans naf_cleannaf et naf_cleannaf_notab*/
29	/* v0.96 M. GASTINEAU 12/01/99 : correction dans le cas ou ICPLX=0 dans */
30	/* naf_modfre et naf_gramsc. */
31	/v0.96 M. GASTINEAU 14/01/99 : ajout du support des fenetres dans /
32	/* naf_mfttab et naf_fftmax. */
33	/v0.97 M. GASTINEAU 26/05/99 : correction bug (0.97/99/05/26/A) dans /
34	/* naf_mftnaf */
35	/* v0.97 M. GASTINEAU 27/05/99 : correction bug (0.97/99/05/27/A) dans le */
36	/* cas ou ICPLX=0 et FS=0 dans naf_modfre et naf_gramsc. */
37
38	#include "complexe.h"
39
40	#define FUNCINLINEEXTERN_MUSTBEINLIB
41
42	/variable globale de naf /
43	t_naf g_NAFVariable;
44
45	/legere difference entre NAF_USE_OPTIMIZE=0 et NAF_USE_OPTIMIZE=1/
46	/* car division differente dans naf_gramsc */
47	#define NAF_USE_OPTIMIZE 1
48
49	/*
50	!-------------------------------------------------------------------------
51	! NAFF.0.84 NUMERICAL ANALYSIS OF FUNDAMENTAL FREQUENCIES
52	! (27 JANVIER 1996)
53	!
54	! (C) JACQUES LASKAR
55	! ASTRONOMIE ET SYSTEMES DYNAMIQUES
56	! BUREAU DES LONGITUDES
57	! 75006 PARIS
58	! EMAIL : LASKAR@BDL.FR
59	!
60	! MAIN REFERENCES :
61	! LASKAR, J.: THE CHAOTIC MOTION OF THE SOLAR SYSTEM. A NUMERICAL
62	! ESTIMATE OF THE SIZE OF THE CHAOTIC ZONES,ICARUS,88,(1990),266--291
63	!
64	! (NAFF.MAC.1.6 - 27/5/98)
65	!
66	!**************************************************************************
67	! THIS PROGRAMM CANNOT BE COPIED,
68	! DISTRIBUTED NOR MODIFIED WITHOUT THE AGREEMENT OF THE AUTHOR.
69	!**************************************************************************
70	!
71	! PROGRAMME D'ANALYSE DE FOURIER AUTOMATIQUE
72	! -- NAFF -- CALCULE UN TERME A LA FOIS ET LE RETRANCHE
73	! A TABS
74	! NOUVELLE VERSION 26/9/87
75	! MODIFIEE POUR NFS LE 18/10/87
76	! MODIFIEE LE 10/4/88 POUR TRAITER
77	! A LA FOIS UN TABLEAU REEL OU COMPLEXE
78	! CALCUL BATCH SUR CIRCE
79	! 8/10/90 MODIF POUR REMETTRE SUR VP (FORTRAN STANDARD)
80	! 13/12/90 MODIF POUR RENDRE PLUS MODULAIRE
81	! 18/12/90 MODIF IL N'Y A PLUS DE DIMENSIONS A REGLER DANS
82	! LES SOUS PROGRAMMES. TABS(2,*) AU LIEU DE TABS(2,0:KTABS)
83	! NFR EST ENCORE EN PARAMETER, MAIS EST RAREMENT A CHANGER.
84	!
85	! IL FAUT REMPLACER LE * DES TABLEAUX
86	! DIMENSIONES T() OU T(N,) PAR T(1) OU T(N,1)
87	! 19/12/90
88	! 16/4/91 COMPILATION SEPAREE DE NAFF EN UNE BIBLIOTHEQUE
89	! 30/4/91 TRAITEMENT DU CAS OU ON RETROUVE LA MEME FREQUENCE
90	! 27/1/96 PLUSIEURS FENETRES
91	! 27/2/97 MODIF POUR DIMENSIONNER TOUT EXPLICITEMENT AVEC
92	! INCLUDE NAFF.INC POUR LES PARAMETRES
93	! 23/45/1997 CORRECTION DE MAXIQUA.F (F. JOUTEL)
94	! 22/4/1998 MODIF POUR QUADRUPLE PRECISION (J. LASKAR)
95	! CHANGEMENT POUR UN EPSILON MACHINE DONNE PAR NAFF.INC
96	! 27/5/98 AMELIORATION DE LA RECHERCHE DU MAXIMUM PAR
97	! LA METHODE DES SECANTES
98	! ROUTINE DE CORRECTION DE LA FREQUENCE PRINCIPALE PAR
99	! DEVELOPPEMENT ASYMTOTIQUE (F. JOUTEL)
100	!***********************************************************************
101	! MODULE NAFF
102	!----------------------------------------------------------------------
103	! PROCEDURE D'UTILISATION:
104	! - PHASE D'INITIALISATION
105	! 1) INITIALISER LES VARIABLES :
106	! DTOUR,KTABS,XH,NTERM,IW,T0,ICPLX,ISEC,NFPRT,IPRT
107	! 2) CALL INITNAF
108	! CALCULE UNIANG,FREFON, EPSM, PI
109	! CREE LES TABLEAUX DYNAMIQUES ((ZTABS)), ((TFS)), ((ZAMP)),
110	! ((ZALP)) et((TWIN))
111	! APPELLE INIWIN
112	! - PHASE D'UTILISATION
113	! 1) EMPLIR LE TABLEAU ZTABS
114	! 2) CALL MFTNAF(NBTERM,EPS)
115	! ANALYSE ((ZTABS)) ET ESSAIE DE RECHERCHER (NBTERM) FREQUENCES,
116	! (EPS) ETANT L'ERREUR TOLEREE
117	! --! MODIFIE ((ZTABS))
118	! - PHASE D'EXPLOITATION OU DE MODIFICATION
119	! * CALL PRTABS(KTABS,ZTABS,IPAS)
120	! IMPRIME LES ELEMENTS DU TABLEAU COMPLEXE (ZTABS(0:KTABS)) TOUS
121	! LES (IPAS), ((IPRT)) DOIT ETRE MIS A 1 AVANT LA PROCEDURE POUR
122	! IMPRIMER QUELQUE CHOSE
123	! * CALL SMOY(ZM)
124	! CALCULE LA MOYENNNE DES ELEMENTS DE ((ZTABS)) ET MODIFIE CE
125	! TABLEAU EN SOUSTRAYANT DE CHACUN DES ELEMENTS CETTE VALEUR
126	! MOYENNE (ZM)
127	! --! MODIFIE ((ZTABS))
128	! * CALL TESSOL (EPS,TFSR,ZAMPR)
129	! TESTE LA COHERENCE DES SOLUTIONS OBTENUES PAR NAFF:
130	! - RECONSTRUIT LE SIGNAL A PARTIR DES ((NFS)) TERMES TROUVES
131	! QUI SONT DANS ((TFS)) ET ((ZAMP))
132	! - ANALYSE CE NOUVEAU SIGNAL, (EPS) ETANT L'ERREUR TOLEREE,
133	! ET RANGE LES NOUVEAUX TERMES DANS (TFSR) ET (ZAMPR)
134	! --! MODIFIE ((ZTABS))
135	! * CALL INIFRE
136	! REMET A ZERO ((TFS)), ((ZAMP)),((ZALP)) et ((NFS))
137	! UTILE QUAND ON BOUCLE SUR PLUSIEURS CAS
138	! * CALL CORRECTION (FREQ)
139	! UTILISE LES TERMES TROUVES POUR CORRIGER LA FREQUENCE
140	! PRINCIPALE, LA NOUVELLE VALEUR CORRIGEE EST (FREQ)
141	! - PHASE DE SORTIE (POUR RECUPERER DE LA PLACE OU POUR REINITIALISER)
142	! 1) CALL CLEANNAF
143	! DESALLOUE LES TABLEAUX DYNAMIQUES ((ZTABS)), ((TFS)), ((ZAMP)),
144	! ((ZALP)) et ((TWIN))
145	!-------------------------------------------------------------------------
146	! VARIABLES INDISPENSABLES
147	! NTERM : LE NOMBRE MAXIMAL DE FREQUENCES A RECHERCHER
148	! KTABS : LE NOMBRE D'INTERVALLES ENTRE LES DONNEES
149	! (IL Y A DONC KTABS+1 DONNEES)
150	! XH : PAS ENTRE DEUX DONNEES
151	! DANS UNE UNITE DE TEMPS QUELCONQUE
152	! LA DUREE DE L'INTERVALLE EST ALORS XH!KTABS
153	! T0 : DATE DE DEPART DES DONNEES
154	! DTOUR : LONGUEUR D'UN TOUR DE CADRAN
155	! ICPLX : 0 SI LA FONCTION EST REELLE, 1 SINON
156	! IW : PARAMETRE DE LA FENETRE D'INTEGRATION
157	! IW=0 PAS DE FENETRE
158	! IW=1 FENETRE DE HANNING etc
159	! ISEC : DRAPEAU DE LA METHODE DES SECANTES
160	! ISEC=1 oui, ISEC=0 non
161	! NFPRT : NUMERO DE LA CONNECTION POUR L'IMPRESSION
162	! IPRT : INDIQUE LE TYPE D'IMPRESSION, EN GENERAL :
163	! -1 PAS D'IMPRESSION
164	! O MESSAGES D'ERREURS EVENTUELLES
165	! 1 RESULTATS
166	! 2 (DEBUG)
167	! TABLEAU OU RANGER LE SIGNAL
168	! ZTABS(0:KTABS) : TABLEAU DE KTABS+1 COMPLEXES
169	!
170	! VARIABLES CALCULEES PAR NAFF
171	! UNIANG : UNITE D'ANGLE = 1RD=UNIANG UNITE D'ANGLE
172	! (SI L'UNITE EST EN SECONDES D'ARC
173	! UNIANG=206264.80624709D0 )
174	! CALCULE EN FONCTION DE DTOUR
175	! FREFON : FREQUENCE FONDAMENTALE
176	! FREFON=2PI/(KTABSXH) OU EN SECONDES
177	! FREFON=360.D03600.D0/(KTABSXH)
178	! NFS : NOMBRE DE FREQUENCES ENTIEREMENT DETERMINEES
179	! AINSI QUE LEUR AMPLITUDE
180	! TFS(NBTERM) : TABLEAU DES FREQUENCES
181	! ZAMP(NBTERM) : AMPLITUDE COMPLEXE DES TERMES
182	! ZALP(NBTERM,NBTERM) : TABLEAU DE CHANGEMENT DE BASE
183	!
184	!-----------------------------------------------------------------------*/
185	static double AF,BF;
186	static double *TWIN=NULL;
187
188	/*!-----------------------------------------------------------------------
189	! VARIABLES A INITIALISER PAR L'UTILISATEUR AVANT DE LANCER INITNAF
190	PUBLIC :: DTOUR,KTABS,XH,NTERM,IW,NFPRT,IPRT;
191	! VARIABLES PUBLIQUES INITIALISEES PAR INITNAF
192	PUBLIC :: EPSM,PI,UNIANG,FREFON;
193	! TABLEAUX PUBLICS CREES PAR INITNAF ET DETRUITES PAR CLEANNAFF
194	PUBLIC :: TFS,ZAMP,ZTABS,ZALP;
195	!-----------------------------------------------------------------------
196	! VARIABLES A INITIALISER PAR L'UTILISATEUR AVANT DE LANCER NAFF
197	PUBLIC :: T0,ICPLX,ISEC
198	! VARIABLES INITIALISES ET UTILISEES PAR NAFF
199	PUBLIC :: NERROR,NFS
200	!-----------------------------------------------------------------------
201	! ROUTINES PUBLIQUES
202	PUBLIC :: INITNAF,CLEANNAF,MFTNAF,PRTABS,SMOY,TESSOL
203	PUBLIC :: INIFRE,CORRECTION*/
204	/v0.96 M. GASTINEAU 06/10/98 : ajout /
205	void naf_initnaf_notab();
206	void naf_cleannaf_notab();
207	/v0.96 M. GASTINEAU 06/10/98 : fin ajout /
208	void naf_initnaf();
209	void naf_inifre();
210	void naf_cleannaf();
211	BOOL naf_mftnaf(int NBTERM, double EPS);
212	void naf_prtabs(int KTABS, t_complexe *ZTABS, int IPAS);
213	void naf_smoy(t_complexe *ZM);
214	BOOL naf_tessol(double EPS, double TFSR, t_complexe ZAMPR);
215	void naf_correction(double *FREQ);
216	void naf_four1(double DATA /tableau commencant a l'indice 1 */,
217	int NN, int ISIGN);
218	void naf_puiss2(int NT, int *N2);
219	/v0.96 M. GASTINEAU 18/12/98 : modification du prototype /
220	/void naf_iniwin();//remplacee par: /
221	void naf_iniwin(double *p_pardTWIN);
222	/v0.96 M. GASTINEAU 18/12/98 : fin modification /
223	void delete_list_fenetre_naf(t_list_fenetre_naf *p_pListFenNaf);
224	t_list_fenetre_naf concat_list_fenetre_naf(t_list_fenetre_naf p_pListFenHead,
225	t_list_fenetre_naf *p_pListFenEnd);
226	t_list_fenetre_naf *cree_list_fenetre_naf(const double p_dFreqMin,
227	const double p_dFreqMax,
228	const int p_iNbTerm);
229
230	/*!-----------------------------------------------------------------------
231	! VARIABLES PRIVEES
232	PRIVATE :: TWIN,AF,BF
233	! ROUTINES PRIVEES
234	PRIVATE :: INIWIN,FRETES,ZTPOW,FFTMAX,FOUR1,PUISS2,MAXX
235	PRIVATE :: MODFRE, GRAMSC,PROSCA,SECANTES,MAXIQUA
236	PRIVATE :: FUNC,FUNCP,PRODER,ZTDER,FREFIN,PROFRE
237	PRIVATE :: ZTPOW2,ZARDYD,PROSCAA,ZTPOW2A,MODTAB*/
238	static void naf_fretes(double FR, int IFLAG, double TOL, int NUMFR);
239	static void naf_ztpow(int N, int N1, t_complexe *ZT, t_complexe ZA, t_complexe ZAST);
240	/* v0.96 M. GASTINEAU 12/01/99 : optimisation */
241	#if NAF_USE_OPTIMIZE==0
242	static void naf_maxx(int N, double T, int INDX);
243	static void naf_fftmax(double *FR);
244	#else /remplacee par:/
245	static int naf_maxx(int N, double *T);
246	static double naf_fftmax(int p_iFrMin, int p_iFrMax, double FREFO2, int KTABS2);
247	#endif /NAF_USE_OPTIMIZE/
248	/* v0.96 M. GASTINEAU 12/01/99 : fin optimisation */
249	static void naf_modtab(int N, double *T) __attribute__((unused));
250	static void naf_modfre(int NUMFR, double A, double B);
251	static BOOL naf_gramsc(double FS, double A, double B);
252	static void naf_prosca(double F1, double F2, t_complexe* ZP)
253	__attribute__((unused));
254	static void naf_proder(double FS, double DER, double A, double B,double RM);
255	static double naf_funcp(double X);
256	static void naf_ztder(int N, int N1, t_complexe ZTF, t_complexe ZTA, double *TW, t_complexe ZA, t_complexe ZAST, double T0, double XH);
257	static void naf_secantes(double X, double PASS, double EPS, double XM, int IPRT, FILE NFPRT);
258	static double naf_func(double X);
259	static void naf_maxiqua(double X, double PASS, double EPS, double XM, double YM, int IPRT, FILE *NFPRT);
260	static void naf_frefin(double FR, double A, double B, double RM, const double RPAS0, const double RPREC);
261	static void naf_ztpow2(int N, int N1, t_complexe ZTF, t_complexe ZTA, double *TW, t_complexe ZA, t_complexe ZAST);
262	static BOOL naf_profre(double FS, double A, double B, double *RMD);
263	static BOOL naf_proscaa(double F1, double F2, t_complexe *ZP);
264	static BOOL naf_zardyd(t_complexe ZT, int N, double H, t_complexe ZOM);
265	static void naf_ztpow2a(int N, int N1, t_complexe ZTF, double TW, t_complexe ZA, t_complexe ZAST);
266
267	/*!------------------------------------------------------------------------
268	CONTAINS*/
269
270
271	/!/
272
273	/*!-----------------------------------------------------------------------
274	!-----------------------------------------------------------------------
275	! ROUTINES DE NAFF
276	!-----------------------------------------------------------------------
277	SUBROUTINE INITNAF
278	IMPLICIT NONE
279	!-----------------------------------------------------------------------
280	! INITNAFF
281	! - effectue les initialisations necessaires : PI,UNIANG,FREFON,EPSM
282	! - alloue les tableaux dynamiques TFS,ZAMP,ZALP,ZTABS,TWIN
283	! - appelle INIWIN
284	!-----------------------------------------------------------------------*/
285	void naf_initnaf(void)
286	{
287	/!----------------- PREMIERES INITIALISATIONS/
288	g_NAFVariable.EPSM = DBL_EPSILON;
289	/PI = ATAN2(1.D0,0.D0)2*/
290	g_NAFVariable.UNIANG = g_NAFVariable.DTOUR/(2*pi) ;
291	g_NAFVariable.FREFON = g_NAFVariable.DTOUR/(g_NAFVariable.KTABS*g_NAFVariable.XH) ;
292	SYSCHECKMALLOCSIZE(g_NAFVariable.TFS, double, g_NAFVariable.NTERM+1);/allocate(TFS(1:NTERM),stat = NERROR)/
293	SYSCHECKMALLOCSIZE(g_NAFVariable.ZAMP, t_complexe, g_NAFVariable.NTERM+1); /allocate(ZAMP(1:NTERM),stat = NERROR)/
294	DIM2(g_NAFVariable.ZALP, (g_NAFVariable.NTERM+1), (g_NAFVariable.NTERM+1), t_complexe,"ZALP"); /allocate(ZALP(1:NTERM,1:NTERM),stat = NERROR)/
295	SYSCHECKMALLOCSIZE(g_NAFVariable.ZTABS, t_complexe, g_NAFVariable.KTABS+1);/* allocate(ZTABS(0:KTABS),stat = NERROR)*/
296	SYSCHECKMALLOCSIZE(TWIN, double, g_NAFVariable.KTABS+1); /allocate(TWIN(0:KTABS),stat = NERROR)/
297	/v0.96 M. GASTINEAU 18/12/98 : modification du prototype /
298	/naf_iniwin(); /
299	naf_iniwin(TWIN);
300	}/* end SUBROUTINE INITNAF*/
301	/*!-----------------------------------------------------------------------
302	subroutine CLEANNAF*/
303	/* v0.96 M. GASTINEAU 06/01/99 : ajout de la liberation de la liste des fenetres */
304	void naf_cleannaf(void)
305	{
306	/*!-----------------------------------------------------------------------
307	! desalloue les tableaux dynamiques
308	!-----------------------------------------------------------------------*/
309
310	SYSFREE(g_NAFVariable.TFS);
311	SYSFREE(g_NAFVariable.ZAMP);
312	HFREE2(g_NAFVariable.ZALP);
313	SYSFREE(g_NAFVariable.ZTABS);
314	SYSFREE(TWIN);
315	/* v0.96 M. GASTINEAU 06/01/99 : ajout */
316	delete_list_fenetre_naf(g_NAFVariable.m_pListFen);
317	g_NAFVariable.m_pListFen =NULL;
318	/* v0.96 M. GASTINEAU 06/01/99 : fin ajout */
319	}/* end subroutine CLEANNAF */
320
321	/-----------------------------------------------------------------------/
322	/* fonction identique a naf_initnaf mais : */
323	/n'alloue pas de tableau pour ZTABS, TFS et ZAMP. /
324	/-----------------------------------------------------------------------/
325	/v0.96 M. GASTINEAU 06/10/98 : ajout /
326	void naf_initnaf_notab()
327	{
328	/!----------------- PREMIERES INITIALISATIONS/
329	g_NAFVariable.EPSM = DBL_EPSILON;
330	g_NAFVariable.UNIANG = g_NAFVariable.DTOUR/(2*pi) ;
331	g_NAFVariable.FREFON = g_NAFVariable.DTOUR/(g_NAFVariable.KTABS*g_NAFVariable.XH) ;
332	DIM2(g_NAFVariable.ZALP, (g_NAFVariable.NTERM+1), (g_NAFVariable.NTERM+1), t_complexe,"ZALP"); /allocate(ZALP(1:NTERM,1:NTERM),stat = NERROR)/
333	SYSCHECKMALLOCSIZE(TWIN, double, g_NAFVariable.KTABS+1); /allocate(TWIN(0:KTABS),stat = NERROR)/
334	/v0.96 M. GASTINEAU 18/12/98 : modification du prototype /
335	/naf_iniwin();/ naf_iniwin(TWIN);
336	}/* end SUBROUTINE naf_initnaf_notab*/
337
338	/-----------------------------------------------------------------------/
339	/* fonction identique a naf_initnaf mais : */
340	/ne libere pas ZTABS, TFS et ZAMP. /
341	/-----------------------------------------------------------------------/
342	/v0.96 M. GASTINEAU 06/10/98 : ajout /
343	/* v0.96 M. GASTINEAU 06/01/99 : ajout de la liberation de la liste des fenetres */
344	void naf_cleannaf_notab()
345	{
346	/*!-----------------------------------------------------------------------
347
348	! desalloue les tableaux dynamiques
349	!-----------------------------------------------------------------------*/
350
351	HFREE2(g_NAFVariable.ZALP);
352	SYSFREE(TWIN);
353	/* v0.96 M. GASTINEAU 06/01/99 : ajout */
354	delete_list_fenetre_naf(g_NAFVariable.m_pListFen);
355	g_NAFVariable.m_pListFen =NULL;
356	/* v0.96 M. GASTINEAU 06/01/99 : fin ajout */
357	}/* end subroutine naf_initnaf_notab */
358
359
360	/*!-----------------------------------------------------------------------
361	SUBROUTINE MFTNAF(NBTERM,EPS)*/
362	/v0.97 M. GASTINEAU 26/05/99 : correction bug 0.97/99/05/26/A pour pas<0 (XH<0) /
363	BOOL naf_mftnaf(int NBTERM, double EPS)
364	{
365	/*!-----------------------------------------------------------------------
366	! MFTNAF
367	! CALCULE UNE APPROXIMATION QUASI PERIODIQUE
368	! DE LA FONCTION TABULEE DANS g_NAFVariable.ZTABS(0:KTABS)
369	!
370	! NBTERM : NOMBRE DE TERMES RECHERCHES (<= NTERM)
371	!
372	!
373	! EPS : PRECISION AVEC LAQUELLE ON RECHERCHE
374	! LES FREQUENCES
375	!
376	!-----------------------------------------------------------------------
377
378	IMPLICIT NONE
379	! (EPS)
380	integer :: NBTERM
381	REAL (8) :: EPS
382	!
383	integer :: I,IFLAG,NUMFR
384	REAL (8) :: TOL,STAREP,FR,A,B,RM
385	!-----------------------------------------------------------------------*/
386	int I, IFLAG, NUMFR = 0;
387	double TOL,STAREP,FR,A,B,RM;
388	/v0.96 M. GASTINEAU 14/01/99 : ajout - support des fenetres/
389	#if NAF_USE_OPTIMIZE>0
390	int iFrMin;
391	int iFrMax;
392	double FREFO2;
393	int KTABS2;
394	#endif /NAF_USE_OPTIMIZE/
395	/v0.96 M. GASTINEAU 14/01/99 : fin ajout/
396
397	if (NBTERM >g_NAFVariable.NTERM)
398	{
399	Myyerror("Nbre de termes cherches trop grand");
400	return FALSE;
401	}
402	TOL = 1.E-4;
403	STAREP=fabs(g_NAFVariable.FREFON)/3;
404	naf_inifre();
405	/v0.96 M. GASTINEAU 14/01/99 : ajout support des fenetres/
406	#if NAF_USE_OPTIMIZE>0
407	naf_puiss2(g_NAFVariable.KTABS+1,&KTABS2);
408	FREFO2=(g_NAFVariable.FREFON*g_NAFVariable.KTABS)/KTABS2;
409	do
410	{
411	if (g_NAFVariable.m_pListFen==NULL)
412	{/pas de fenetre => on recherche dans toutes les frequences /
413	iFrMin=KTABS2;
414	iFrMax=KTABS2;
415	}
416	else
417	{
418	const int iMaxValue=KTABS2/2-1;
419	NBTERM=g_NAFVariable.m_pListFen->iNbTerme;
420	iFrMin=(int)(g_NAFVariable.m_pListFen->dFreqMin/FREFO2);
421	iFrMax=(int)(g_NAFVariable.m_pListFen->dFreqMax/FREFO2);
422	/v0.97 M. GASTINEAU 26/05/99 : ajout - correction bug (cas pas <0) /
423	if (iFrMin>iFrMax)
424	{/swap/
425	double temp = iFrMin;
426	iFrMin = iFrMax;
427	iFrMax = (int)temp;
428	}
429	/v0.97 M. GASTINEAU 26/05/99 : fin ajout /
430	if (iFrMin<-iMaxValue)
431	{ iFrMin=-iMaxValue; }
432	else if (iFrMin>iMaxValue)
433	{ iFrMin=+iMaxValue; }
434	if (iFrMax<-iMaxValue)
435	{ iFrMax=-iMaxValue; }
436	else if (iFrMax>iMaxValue)
437	{ iFrMax=+iMaxValue; }
438	}
439	#endif /NAF_USE_OPTIMIZE/
440	/v0.96 M. GASTINEAU 14/01/99 : fin ajout /
441	for(I=1;I<=NBTERM; I++)
442	{
443	/v0.96 M. GASTINEAU 14/01/99 : modification support des fenetres/
444	#if NAF_USE_OPTIMIZE==0
445	naf_fftmax(&FR);
446	#else /remplacee par:/
447	FR=naf_fftmax(iFrMin,iFrMax,FREFO2,KTABS2);
448	#endif /NAF_USE_OPTIMIZE/
449	/v0.96 M. GASTINEAU 14/01/99 : fin modification/
450	naf_frefin(&FR,&A,&B,&RM,STAREP,EPS);
451	naf_fretes(FR,&IFLAG,TOL,&NUMFR);
452	if (IFLAG == 0) break; /GOTO 999/
453	if (IFLAG == 1)
454	{
455	if(naf_gramsc(FR,A,B)==FALSE)
456	{
457	return FALSE;
458	}
459	g_NAFVariable.NFS++; /g_NAFVariable.NFS=g_NAFVariable.NFS+1;/
460	}
461	if (IFLAG==-1)
462	{
463	naf_modfre(NUMFR,&A,&B);
464	}
465	}
466	/v0.96 M. GASTINEAU 14/01/99 : ajout - support des fenetres/
467	#if NAF_USE_OPTIMIZE>0
468	/passage a la fenetre suivante /
469	if (g_NAFVariable.m_pListFen!=NULL)
470	{
471	t_list_fenetre_naf *pListTemp=g_NAFVariable.m_pListFen;
472	g_NAFVariable.m_pListFen=pListTemp->suivant;
473	SYSFREE(pListTemp);
474	}
475	} while (g_NAFVariable.m_pListFen!=NULL);
476	#endif /NAF_USE_OPTIMIZE/
477	/v0.96 M. GASTINEAU 14/01/99 : fin ajout/
478
479	/999 CONTINUE/
480	return TRUE;
481	} /* END SUBROUTINE MFTNAF */
482	/*!
483	SUBROUTINE PRTABS(KTABS,g_NAFVariable.ZTABS,IPAS)*/
484	void naf_prtabs(int KTABS, t_complexe *ZTABS, int IPAS)
485	{
486	/*
487	!-----------------------------------------------------------------------
488	! IMPRESSION DE g_NAFVariable.ZTABS
489	!-----------------------------------------------------------------------
490
491	IMPLICIT NONE
492	! (KTABS,g_NAFVariable.ZTABS,IPAS)
493	integer :: KTABS,IPAS
494	complex (8) :: g_NAFVariable.ZTABS(0:KTABS)
495	!
496	integer :: I
497	!*/
498	int I;
499
500	if (g_NAFVariable.IPRT==1)
501	{
502	for (I=0;I<=KTABS;I+=IPAS)
503	{
504	fprintf(g_NAFVariable.NFPRT,"%6d %-+20.15E %-+20.15E\n",I, ZTABS[I].reel, ZTABS[I].imag);
505	/* WRITE(g_NAFVariable.NFPRT,1000) I,DREAL(ZTABS(I)),DIMAG(ZTABS(I))*/
506	}
507	}
508	/1000 FORMAT (1X,I6,2X,2D20.6)/
509	}/* END SUBROUTINE PRTABS */
510	/*!
511	!
512	SUBROUTINE SMOY(ZM)*/
513	void naf_smoy(t_complexe *ZM)
514	{
515	/*
516	!-----------------------------------------------------------------------
517	! CALCUL DES MOYENNNES DE g_NAFVariable.ZTABS ET SOUSTRAIT DE CHACUNE
518	! LA VALEUR MOYENNE
519
520	!-----------------------------------------------------------------------
521
522	IMPLICIT NONE
523	! (ZM)
524	complex (8) :: ZM
525	!
526	integer :: I
527	!-----------------------------------------------------------------------
528	*/
529	int I;
530	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
531	#if NAF_USE_OPTIMIZE==0
532	*ZM=cmplx(0.E0,0.E0);
533	for (I=0; I<=g_NAFVariable.KTABS; I++)
534	{
535	ZM = addcomplexe(ZM,g_NAFVariable.ZTABS[I]);
536	}
537	ZM=muldoublcomplexe(1E0/((double)(g_NAFVariable.KTABS+1)), ZM); /ZM=ZM/(g_NAFVariable.KTABS+1)/
538	for (I=0; I<=g_NAFVariable.KTABS; I++)
539	{
540	g_NAFVariable.ZTABS[I]=subcomplexe(g_NAFVariable.ZTABS[I],*ZM);
541	}
542	#else /remplacees par: /
543	t_complexe *pzarTabs;
544	double dNbKTabs1=g_NAFVariable.KTABS+1;
545	i_compl_cmplx(ZM,0.E0,0.E0);
546	for (I=0, pzarTabs = g_NAFVariable.ZTABS;
547	I<=g_NAFVariable.KTABS;
548	I++, pzarTabs++)
549	{
550	i_compl_padd(ZM,pzarTabs);
551	}
552	i_compl_pdivdoubl(ZM,&dNbKTabs1);/ZM=ZM/(g_NAFVariable.KTABS+1)/
553	for (I=0, pzarTabs = g_NAFVariable.ZTABS;
554	I<=g_NAFVariable.KTABS;
555	I++, pzarTabs++)
556	{
557	i_compl_psub(pzarTabs,ZM);
558	}
559	#endif /NAF_USE_OPTIMIZE/
560	/* v0.96 M. GASTINEAU 01/12/98 : fin modification */
561	}/* END SUBROUTINE SMOY*/
562	/*!
563	!
564	SUBROUTINE FRETES (FR,IFLAG,TOL,NUMFR)*/
565	void naf_fretes(double FR, int IFLAG, double TOL, int NUMFR)
566	{
567	/!*********************************************************************
568	! TEST DE LA NOUVELLE FREQUENCE TROUVEE PAR RAPPORT AUX ANCIENNES.
569	! LA DISTANCE ENTRE DEUX FREQ DOIT ETRE DE g_NAFVariable.FREFON
570	!
571	! RENVOIE IFLAG = 1 SI LE TEST REUSSIT (ON PEUT CONTINUER)
572	! IFLAG = 0 SI LE TEST ECHOUE (IL VAUT MIEUX S'ARRETER)
573	! IFLAG = -1 SI TEST < ECART, MAIS TEST/ECART < TOL
574	! (ON RETROUVE PRATIQUEMENT LA MEME FREQUENCE
575	! D'INDICE NFR)
576	! TOL (ENTREE) TOLERANCE ( 1.D-7) EST UN BON EXEMPLE
577	! NFR (SORTIE) INDICE DE LA FREQUENCE RETROUVEE
578	!
579	!***********************************************************************
580
581	IMPLICIT NONE
582	! (FR,IFLAG,TOL,NUMFR)
583	integer :: IFLAG,NUMFR
584	REAL (8) :: FR,TOL
585	!
586	integer :: I
587	REAL (8) :: ECART,TEST
588	! */
589	int I;
590	double ECART, TEST;
591	*IFLAG = 1;
592	ECART = fabs(g_NAFVariable.FREFON) ;
593	for (I = 1; I<=g_NAFVariable.NFS; I++)
594	{
595	TEST = fabs(g_NAFVariable.TFS[I] - FR);
596	if (TEST<ECART)
597	{
598	if ((TEST/ECART)<TOL)
599	{
600	*IFLAG = -1;
601	*NUMFR = I;
602	if (g_NAFVariable.IPRT>=1)
603	{
604	fprintf(g_NAFVariable.NFPRT, "TEST/ECART = %g ON CONTINUE\n", TEST/ECART);
605	}
606	break; /GOTO 999/
607	}
608	else
609	{
610
611	*IFLAG = 0 ;
612	if (g_NAFVariable.IPRT>=0)
613	{
614	fprintf(g_NAFVariable.NFPRT,"TEST = %g ECART = %g \n", TEST, ECART);
615	fprintf(g_NAFVariable.NFPRT,"FREQUENCE FR = %g TROP PROCHE DE %g\n", FR, g_NAFVariable.TFS[I]);
616	}
617	}
618	}
619	}
620	/999 CONTINUE/
621	}/* END SUBROUTINE FRETES*/
622	/*!
623	SUBROUTINE ZTPOW (N,N1,ZT,ZA,ZAST)*/
624	void naf_ztpow(int N, int N1, t_complexe *ZT, t_complexe ZA, t_complexe ZAST)
625	{
626	/*!-----------------------------------------------------------------------
627	! ZTPOW CALCULE ZT(I) = ZASTZA*I EN VECTORIEL
628	! ZT(N)
629	!-----------------------------------------------------------------------
630	IMPLICIT NONE
631	! (N,N1,ZT,ZA,ZAST)
632	integer :: N,N1
633	complex (8) :: ZT(0:N),ZA,ZAST
634	!
635	integer :: I,INC,NX,IT,NT
636	complex (8) ZT1,ZINC
637	! */
638	int I,INC,NX,IT,NT;
639	t_complexe ZT1,ZINC;
640	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
641	#if NAF_USE_OPTIMIZE>0
642	t_complexe *pzarZT;
643	#endif /NAF_USE_OPTIMIZE/
644
645	if (N<N1-1)
646	{
647	printf("DANS ZTPOW, N = %d\n", N);
648	return;
649	}
650	/!----------! /
651	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
652	#if NAF_USE_OPTIMIZE==0
653	ZT[0] = mulcomplexe(ZAST,ZA);
654	for(I = 1; I<N1; I++)
655	{
656	ZT[I] = mulcomplexe(ZT[I-1],ZA);
657	}
658	ZT1=divcomplexe(ZT[N1-1], ZAST); /ZT1 = ZT(N1-1)/ZAST/
659	ZINC=cmplx(1E0, 0E0); /ZINC= 1;/
660	INC =0;
661	NT = (N+1)/N1;
662	for (IT = 2; IT<= NT; IT++)
663	{
664	ZINC = mulcomplexe(ZINC,ZT1);
665	INC += N1; /INC = INC + N1/
666	for (I = 0; I<N1; I++)
667	{
668	ZT[INC+I]=mulcomplexe(ZT[I], ZINC); /ZT(INC +I) = ZT(I)ZINC*/
669	}
670	}
671	ZINC = mulcomplexe(ZINC, ZT1); /ZINC = ZINCZT1*/
672	INC += N1; /INC = INC + N1;/
673	NX = N+1-NT*N1;
674	for (I = 0; I<NX; I++)
675	{
676	ZT[INC+I]=mulcomplexe(ZT[I], ZINC); /ZT(INC +I) = ZT(I)ZINC*/
677	}
678	#else /remplacee par: /
679	ZT[0] = i_compl_mul(ZAST,ZA);
680	for(I = 1, pzarZT=ZT+1 ; I<N1; I++,pzarZT++)
681	{
682	pzarZT = i_compl_mul((pzarZT-1),ZA);
683	}
684	ZT1=i_compl_div(ZT[N1-1], ZAST); /ZT1 = ZT(N1-1)/ZAST/
685	i_compl_cmplx(&ZINC,1E0, 0E0); /ZINC= 1;/
686	INC =0;
687	NT = (N+1)/N1;
688	for (IT = 2; IT<= NT; IT++)
689	{
690	i_compl_pmul(&ZINC,&ZT1);
691	INC += N1; /INC = INC + N1/
692	for (I = 0; I<N1; I++)
693	{
694	ZT[INC+I]=i_compl_mul(ZT[I], ZINC); /ZT(INC +I) = ZT(I)ZINC*/
695	}
696	}
697	i_compl_pmul(&ZINC, &ZT1); /ZINC = ZINCZT1*/
698	INC += N1; /INC = INC + N1;/
699	NX = N+1-NT*N1;
700	for (I = 0; I<NX; I++)
701	{
702	ZT[INC+I]=i_compl_mul(ZT[I], ZINC); /ZT(INC +I) = ZT(I)ZINC*/
703	}
704	#endif /NAF_USE_OPTIMIZE/
705	/* v0.96 M. GASTINEAU 01/12/98 : fin modification */
706	}/* END SUBROUTINE ZTPOW*/
707	/*!
708	SUBROUTINE TESSOL (EPS,TFSR,ZAMPR)*/
709	BOOL naf_tessol(double EPS, double TFSR, t_complexe ZAMPR)
710	{
711	/*!-----------------------------------------------------------------------
712	! TESSOL
713	! SOUS PROGRAMME POUR VERIFIER LA PRECISION DES
714	! SOLUTIONS OBTENUES PAR ANALYSE DE FOURIER
715	! ON ANALYSE A NOUVEAU LA SOLUTION ET ON COMPARE
716	! LES TERMES
717	!
718	! J. LASKAR 25/5/89
719	!-----------------------------------------------------------------------
720
721	IMPLICIT NONE
722	! (EPS,TFSR,ZAMPR)
723	REAL (8) :: EPS, TFSR(g_NAFVariable.NTERM)
724	complex (8) :: ZAMPR(g_NAFVariable.NTERM)
725	!
726	integer :: IT,IFR,JFR, NVTERM
727	REAL (8) :: OFFSET
728	complex (8) :: ZI,ZA,ZOM,ZEX
729	REAL (8), dimension(:),allocatable :: TFST
730	complex (8),allocatable,dimension(:) :: ZAMPT
731	complex (8),allocatable,dimension(:,:) :: ZALPT
732	complex (8) :: ZINC
733	complex (8), dimension (:), allocatable :: ZT
734	*/
735	int IT,IFR,JFR, NVTERM;
736	double OFFSET;
737	t_complexe ZI,ZA,ZOM,ZEX;
738	double *TFST=NULL;
739	t_complexe *ZAMPT=NULL;
740	t_complexe **ZALPT=NULL;
741	t_complexe ZINC;
742	t_complexe *ZT=NULL;
743	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
744	#if NAF_USE_OPTIMIZE>0
745	t_complexe *pzarTab;
746	#endif /**/
747	SYSCHECKMALLOCSIZE(ZAMPT, t_complexe, g_NAFVariable.NTERM+1); /* allocate(ZAMPT(1:g_NAFVariable.NTERM),stat = NERROR)*/
748	DIM2(ZALPT, (g_NAFVariable.NTERM+1), (g_NAFVariable.NTERM+1), t_complexe, "ZALPT"); /allocate(ZALPT(1:g_NAFVariable.NTERM,1:g_NAFVariable.NTERM),stat = NERROR)/
749	SYSCHECKMALLOCSIZE(TFST, double, g_NAFVariable.NTERM+1);/* allocate(TFST(1:g_NAFVariable.NTERM),stat = NERROR)*/
750	SYSCHECKMALLOCSIZE(ZT, t_complexe, g_NAFVariable.KTABS+1);/* allocate (ZT (0:g_NAFVariable.KTABS),stat = NERROR)*/
751
752	/!************************************************************/
753	OFFSET=g_NAFVariable.KTABS*g_NAFVariable.XH;
754	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
755	#if NAF_USE_OPTIMIZE==0
756	/!----------! INITIALISATION DE g_NAFVariable.ZTABS/
757	for ( IT=0;IT<=g_NAFVariable.KTABS; IT++)
758	{
759	g_NAFVariable.ZTABS[IT]=cmplx(0.E0,0.E0);
760	}
761	/!----------! CALCUL DE LA NOUVELLE SOLUTION/
762	ZI=cmplx(0.E0,1.E0);
763	for ( IFR = 1; IFR<=g_NAFVariable.NFS; IFR++)
764	{
765	ZOM = muldoublcomplexe(g_NAFVariable.TFS[IFR]/g_NAFVariable.UNIANG,ZI); /* ZOM=g_NAFVariable.TFS(IFR)/g_NAFVariable.UNIANGZI /
766	ZA=g_NAFVariable.ZAMP[IFR];
767	/*!-----------! ATTENTION ICI (CAS REEL) ON AURA AUSSI LE TERME CONJUGUE
768	!-----------! QU'ON NE CALCULE PAS. LE TERME TOTAL EST
769	!-----------! 2RE(g_NAFVariable.ZAMP(I)EXP(ZIg_NAFVariable.TFS(I)T) )
770	!-----------! ON RAJOUTE LA CONTRIBUTION TOTALE DU TERME g_NAFVariable.TFS(I) DANS TAB2*/
771	if (g_NAFVariable.ICPLX==1)
772	{
773	ZEX = mulcomplexe(ZA,expcomplexe(muldoublcomplexe(g_NAFVariable.T0+OFFSET-g_NAFVariable.XH,ZOM))); /ZEX = ZAEXP ((g_NAFVariable.T0+OFFSET-g_NAFVariable.XH)ZOM)/
774	ZINC= expcomplexe(muldoublcomplexe(g_NAFVariable.XH,ZOM)); /ZINC= EXP (g_NAFVariable.XHZOM)*/
775	naf_ztpow(g_NAFVariable.KTABS,64,ZT,ZINC,ZEX); /naf_ztpow(g_NAFVariable.KTABS+1,64,ZT,ZINC,ZEX);/
776	for (IT=0; IT<=g_NAFVariable.KTABS; IT++)
777	{
778	g_NAFVariable.ZTABS[IT]=addcomplexe(g_NAFVariable.ZTABS[IT], ZT[IT]); /g_NAFVariable.ZTABS(IT)=g_NAFVariable.ZTABS(IT)+ZT(IT)/
779	}
780	}
781	else
782	{
783	ZEX = mulcomplexe(ZA,expcomplexe(muldoublcomplexe(g_NAFVariable.T0+OFFSET-g_NAFVariable.XH,ZOM))); /ZEX = ZAEXP ((g_NAFVariable.T0+OFFSET-g_NAFVariable.XH)ZOM)/
784	ZINC= expcomplexe(muldoublcomplexe(g_NAFVariable.XH,ZOM)); /ZINC= EXP (g_NAFVariable.XHZOM)*/
785	naf_ztpow(g_NAFVariable.KTABS,64,ZT,ZINC,ZEX); /naf_ztpow(g_NAFVariable.KTABS+1,64,ZT,ZINC,ZEX);/
786	for (IT=0; IT<=g_NAFVariable.KTABS; IT++)
787	{
788	g_NAFVariable.ZTABS[IT].reel += ZT[IT].reel; /g_NAFVariable.ZTABS(IT)=g_NAFVariable.ZTABS(IT)+DREAL(ZT(IT))/
789	}
790	}
791	}
792	#else /remplacee par:/
793	/!----------! INITIALISATION DE g_NAFVariable.ZTABS/
794	for ( IT=0, pzarTab=g_NAFVariable.ZTABS;
795	IT<=g_NAFVariable.KTABS;
796	IT++, pzarTab++)
797	{
798	i_compl_cmplx(pzarTab, 0.E0, 0.E0);
799	}
800	/!----------! CALCUL DE LA NOUVELLE SOLUTION/
801	i_compl_cmplx(&ZI,0.E0,1.E0);
802	for ( IFR = 1; IFR<=g_NAFVariable.NFS; IFR++)
803	{
804	ZOM = i_compl_muldoubl(g_NAFVariable.TFS[IFR]/g_NAFVariable.UNIANG,ZI); /* ZOM=g_NAFVariable.TFS(IFR)/g_NAFVariable.UNIANGZI /
805	ZA=g_NAFVariable.ZAMP[IFR];
806	/*!-----------! ATTENTION ICI (CAS REEL) ON AURA AUSSI LE TERME CONJUGUE
807	!-----------! QU'ON NE CALCULE PAS. LE TERME TOTAL EST
808	!-----------! 2RE(g_NAFVariable.ZAMP(I)EXP(ZIg_NAFVariable.TFS(I)T) )
809	!-----------! ON RAJOUTE LA CONTRIBUTION TOTALE DU TERME g_NAFVariable.TFS(I) DANS TAB2*/
810	if (g_NAFVariable.ICPLX==1)
811	{
812	ZEX = i_compl_mul(ZA,i_compl_exp(i_compl_muldoubl(g_NAFVariable.T0+OFFSET-g_NAFVariable.XH,ZOM))); /ZEX = ZAEXP ((g_NAFVariable.T0+OFFSET-g_NAFVariable.XH)ZOM)/
813	ZINC= i_compl_exp(i_compl_muldoubl(g_NAFVariable.XH,ZOM)); /ZINC= EXP (g_NAFVariable.XHZOM)*/
814	naf_ztpow(g_NAFVariable.KTABS,64,ZT,ZINC,ZEX); /naf_ztpow(g_NAFVariable.KTABS+1,64,ZT,ZINC,ZEX);/
815	for (IT=0, pzarTab=g_NAFVariable.ZTABS;
816	IT<=g_NAFVariable.KTABS;
817	IT++, pzarTab++)
818	{
819	i_compl_padd(pzarTab,ZT+IT); /g_NAFVariable.ZTABS(IT)=g_NAFVariable.ZTABS(IT)+ZT(IT)/
820	}
821	}
822	else
823	{
824	ZEX = i_compl_mul(ZA,i_compl_exp(i_compl_muldoubl(g_NAFVariable.T0+OFFSET-g_NAFVariable.XH,ZOM))); /ZEX = ZAEXP ((g_NAFVariable.T0+OFFSET-g_NAFVariable.XH)ZOM)/
825	ZINC= i_compl_exp(i_compl_muldoubl(g_NAFVariable.XH,ZOM)); /ZINC= EXP (g_NAFVariable.XHZOM)*/
826	naf_ztpow(g_NAFVariable.KTABS,64,ZT,ZINC,ZEX); /naf_ztpow(g_NAFVariable.KTABS+1,64,ZT,ZINC,ZEX);/
827	for (IT=0; IT<=g_NAFVariable.KTABS; IT++)
828	{
829	g_NAFVariable.ZTABS[IT].reel += ZT[IT].reel; /g_NAFVariable.ZTABS(IT)=g_NAFVariable.ZTABS(IT)+DREAL(ZT(IT))/
830	}
831	}
832	}
833	#endif /NAF_USE_OPTIMIZE==0/
834	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
835	/!----------------------/
836	SYSFREE(ZT); /deallocate(ZT)/
837	/!-----------! SAUVEGARDE DANS DES TABLEAUX ANNEXES DE LA SOLUTION/
838	for( IFR = 1; IFR<=g_NAFVariable.NFS; IFR++)
839	{
840	TFST[IFR] = g_NAFVariable.TFS[IFR];
841	ZAMPT[IFR] = g_NAFVariable.ZAMP[IFR];
842	for(JFR = 1;JFR <=g_NAFVariable.NFS; JFR++)
843	{
844	ZALPT[IFR][JFR] = g_NAFVariable.ZALP[IFR][JFR];
845	}
846	}
847	/!-----------! NOUVEAU CALCUL DE LA SOLUTION/
848	NVTERM=g_NAFVariable.NFS;
849	/!*DEBUG CALL PRTABS(g_NAFVariable.KTABS,g_NAFVariable.ZTABS,g_NAFVariable.KTABS/10)/
850	if (naf_mftnaf (NVTERM,EPS)==FALSE)
851	{
852	HFREE2(ZALPT);
853	SYSFREE(ZAMPT);
854	SYSFREE(TFST);
855	return FALSE;
856	}
857	/!-----------! RESULTATS/
858	for( IFR = 1; IFR<=g_NAFVariable.NFS; IFR++)
859	{
860	TFSR[IFR] = g_NAFVariable.TFS[IFR];
861	ZAMPR[IFR] = g_NAFVariable.ZAMP[IFR];
862	}
863	/!-----------! RESTITUTION DE LA SOLUTION/
864	for( IFR = 1; IFR<=g_NAFVariable.NFS; IFR++)
865	{
866	g_NAFVariable.TFS[IFR] = TFST[IFR];
867	g_NAFVariable.ZAMP[IFR] = ZAMPT[IFR];
868	for(JFR = 1;JFR <=g_NAFVariable.NFS; JFR++)
869	{
870	g_NAFVariable.ZALP[IFR][JFR] = ZALPT[IFR][JFR];
871	}
872	}
873	HFREE2(ZALPT);
874	SYSFREE(ZAMPT);
875	SYSFREE(TFST);
876	return TRUE;
877	}/* END SUBROUTINE TESSOL*/
878	/*!
879	!
880	SUBROUTINE INIFRE*/
881	void naf_inifre()
882	{
883	/*
884	!************************************************************************
885	! REMET A ZERO g_NAFVariable.TFS, g_NAFVariable.ZAMP,g_NAFVariable.ZALP et g_NAFVariable.NFS
886	! UTILE QUAND ON BOUCLE SUR PLUSIEURS CAS
887	!
888	!***********************************************************************
889
890	IMPLICIT NONE
891	*/
892	int I,J;
893	/* v0.96 M. GASTINEAU 12/01/99 : optimisation */
894	#if NAF_USE_OPTIMIZE==0
895	t_complexe cZERO;
896	cZERO=cmplx(0.E0,0.E0);
897	for(I = 1;I<= g_NAFVariable.NTERM; I++)
898	{
899	g_NAFVariable.TFS[I] = 0.E0;
900	g_NAFVariable.ZAMP[I] = cZERO;
901	for(J = 1; J<= g_NAFVariable.NTERM; J++)
902	{
903	g_NAFVariable.ZALP[I][J] = cZERO;
904	}
905	}
906	#else /* remplace par: */
907	const int iNterm=g_NAFVariable.NTERM;
908	const double dZero=0.E0;
909	double *pdTFS=g_NAFVariable.TFS+1;
910	double pdZAMP=(double)(g_NAFVariable.ZAMP+1);
911	double *pdZALP;
912	for(I = 1;I<= iNterm; I++)
913	{
914	pdZALP=(double*)(g_NAFVariable.ZALP[I]);
915	*pdTFS++ = dZero;
916	*pdZAMP++ = dZero;
917	*pdZAMP++ = dZero;
918	for(J = 1; J<= iNterm; J++)
919	{
920	*pdZALP++ = dZero;
921	*pdZALP++ = dZero;
922	}
923	}
924
925	#endif /NAF_USE_OPTIMIZE/
926	/* v0.96 M. GASTINEAU 12/01/99 : fin optimisation */
927	g_NAFVariable.NFS = 0;
928	} /* END SUBROUTINE INIFRE*/
929	/*!
930	!
931	SUBROUTINE FFTMAX(FR)*/
932	/*
933	!-----------------------------------------------------------------------
934	! FFTMAX MODIF NICE 90 (J. LASKAR)
935	! MODIF 13/12/90
936	! VERSION STANDARD. (FOUR1)
937	! ANALYSE DE FOURIER RAPIDE
938	! RECHERCHE D'UN MAXIMUM DE SPECTRE POUR UNE FREQ. F
939	!***********************************************************************
940	! ON UTILISERA CELA ULTERIEUREMENT (13/12/90)
941	! NFR11, NFR22 NUMERO MIN ET MAX DES FREQUENCES RECHERCHEES
942	! NFR1 ET NFR2 INCLUSES
943	! -NRTAB<NFR1<NFR2<NRTAB
944	!
945	! RTAB(NRTAB,2) EN 1 AMPLITUDE DE SFREQUENCES POSITIVES
946	! 2 AMPLITUDE DES FREQUENCES NEGATIVES
947	! DISTANCE ENTRE DEUX LIGNES g_NAFVariable.FREFON
948	!-----------------------------------------------------------------------
949
950	IMPLICIT NONE
951	! (FR)
952	REAL (8) :: FR
953	!
954	integer :: KTABS2,ISG,IPAS,I,IDV,INDX,IFR
955	REAL (8) :: FREFO2
956	REAL (8), dimension(:),allocatable :: TAB
957	REAL (8), dimension(:),allocatable :: RTAB
958	! */
959	/* v0.96 M. GASTINEAU 14/01/99 : support des fenetres */
960	#if NAF_USE_OPTIMIZE==0
961	void naf_fftmax(double *FR)
962	{
963	int KTABS2,ISG,IPAS,I,IDV,INDX,IFR;
964	double FREFO2;
965	double pdTAB=NULL; /=TAB*/
966	double *RTAB=NULL;
967
968	naf_puiss2(g_NAFVariable.KTABS+1,&KTABS2);
969	/! /
970	SYSCHECKMALLOCSIZE(pdTAB,double,2KTABS2); pdTAB--;/ allocate(TAB(2KTABS2),stat = NERROR)/
971	SYSCHECKMALLOCSIZE(RTAB,double,KTABS2);/allocate(RTAB(0:KTABS2-1),stat = NERROR)/
972	/!/
973	FREFO2=(g_NAFVariable.FREFON*g_NAFVariable.KTABS)/KTABS2;
974	/!***************** */
975	if (g_NAFVariable.IPRT==1)
976	{
977	fprintf(g_NAFVariable.NFPRT,"KTABS2= %d FREFO2= %g\n",KTABS2, FREFO2);
978	}
979	/!**************! CALCUL DES FREQUENCES /
980	ISG=-1;
981	IDV=KTABS2;
982	IPAS=1;
983	for(I=0;I<=KTABS2-1; I++)
984	{
985	pdTAB[2I+1]=g_NAFVariable.ZTABS[I].reelTWIN[I];/pdTAB(2I+1)=DREAL(g_NAFVariable.ZTABS(I))TWIN(I)/
986	pdTAB[2I+2]=g_NAFVariable.ZTABS[I].imagTWIN[I]; /pdTAB(2I+2)=DIMAG(g_NAFVariable.ZTABS(I))TWIN(I)/
987	}
988	naf_four1(pdTAB,KTABS2,ISG);
989	for(I=0;I<=KTABS2-1; I++)
990	{
991	RTAB[I]=sqrt(pdTAB[2I+1]pdTAB[2I+1]+pdTAB[2I+2]pdTAB[2I+2])/IDV;
992	}
993	pdTAB++; SYSFREE(pdTAB);
994	/!********************* CALL MODTAB(KTABS2,RTAB)*/
995	naf_maxx(KTABS2-1, RTAB, &INDX);/naf_maxx(KTABS2,RTAB,&INDX);/
996	/* car naf_maxx travaille sur des tableaux de 1 a N */
997
998	IFR=((INDX+1)<=(KTABS2/2))?INDX:INDX-KTABS2; /*IF (INDX.LE.KTABS2/2) THEN
999	IFR = INDX-1
1000	ELSE
1001	IFR = INDX-1-KTABS2
1002	ENDIF*/
1003
1004	FR = IFRFREFO2*IPAS;
1005	if (g_NAFVariable.IPRT==1)
1006	{
1007	fprintf(g_NAFVariable.NFPRT,"IFR=%d FR=%g RTAB=%g INDX=%d KTABS2=%d\n",IFR,*FR, RTAB[INDX],INDX,KTABS2);
1008	}
1009	SYSFREE(RTAB);
1010	}
1011	#else /remplacee par:/
1012	double naf_fftmax(int p_iFrMin, int p_iFrMax, double FREFO2, int KTABS2)
1013	{
1014	/* la fonction retourne la frequence FR dertminee */
1015	/* On suppose p_iFrMin < p_iFrMax */
1016	double FR;
1017	int ISG,IPAS,I,INDX,IFR;
1018	/double FREFO2;/
1019	double pdTAB=NULL; /=TAB*/
1020	double *RTAB=NULL;
1021	int iKTABS2m1, iKTABS2;
1022	double dDIV;
1023	double pdTABTemp1,pdTABTemp2;
1024
1025
1026
1027	/naf_puiss2(g_NAFVariable.KTABS+1,&KTABS2);//v0.96 M. GASTINEAU 14/01/99 /
1028	iKTABS2 = KTABS2;
1029	iKTABS2m1 = iKTABS2-1;
1030	/! /
1031	SYSCHECKMALLOCSIZE(pdTAB,double,2iKTABS2); / allocate(TAB(2KTABS2),stat = NERROR)/
1032	SYSCHECKMALLOCSIZE(RTAB,double,iKTABS2);/allocate(RTAB(0:KTABS2-1),stat = NERROR)/
1033	/!/
1034	/* FREFO2=(g_NAFVariable.FREFONg_NAFVariable.KTABS)/iKTABS2;//v0.96 M. GASTINEAU 14/01/99 /
1035	/!***************** */
1036	if (g_NAFVariable.IPRT==1)
1037	{
1038	fprintf(g_NAFVariable.NFPRT,"KTABS2= %d FREFO2= %g\n",iKTABS2, FREFO2);
1039	}
1040	/!**************! CALCUL DES FREQUENCES /
1041	ISG=-1;
1042	dDIV=iKTABS2;
1043	IPAS=1;
1044	for(I=0, pdTABTemp1=pdTAB, pdTABTemp2=(double*)(g_NAFVariable.ZTABS);
1045	I<=iKTABS2m1;
1046	I++)
1047	{
1048	pdTABTemp1++ = (pdTABTemp2++) * TWIN[I];/pdTAB(2I+1)=DREAL(g_NAFVariable.ZTABS(I))TWIN(I)/
1049	pdTABTemp1++ = (pdTABTemp2++) * TWIN[I]; /pdTAB(2I+2)=DIMAG(g_NAFVariable.ZTABS(I))TWIN(I)/
1050	}
1051	naf_four1(pdTAB-1,iKTABS2,ISG);
1052	for(I=0, pdTABTemp1=pdTAB, pdTABTemp2=pdTABTemp1+1;
1053	I<=iKTABS2m1;
1054	I++, pdTABTemp1+=2, pdTABTemp2+=2)
1055	{
1056	RTAB[I]=sqrt((pdTABTemp1)(pdTABTemp1)+(pdTABTemp2)(pdTABTemp2))/dDIV;
1057	}
1058	SYSFREE(pdTAB);
1059	/!********************* CALL MODTAB(KTABS2,RTAB)*/
1060	/v0.96 M. GASTINEAU 14/01/99 : modification pour le support des fenetres /
1061	/INDX=naf_maxx(iKTABS2m1, RTAB);//naf_maxx(KTABS2,RTAB,&INDX);/
1062	/IFR=((INDX+1)<=(iKTABS2/2))?INDX:INDX-iKTABS2;/ /*IF (INDX.LE.KTABS2/2) THEN
1063	IFR = INDX-1
1064	ELSE
1065	IFR = INDX-1-KTABS2
1066	ENDIF*/
1067	/* FR = IFRFREFO2*IPAS;
1068	if (g_NAFVariable.IPRT==1)
1069	{
1070	fprintf(g_NAFVariable.NFPRT,"IFR=%d FR=%g RTAB=%g INDX=%d KTABS2=%d\n",IFR,*FR, RTAB[INDX],INDX,iKTABS2);
1071	}
1072	SYSFREE(RTAB); */
1073	/remplacee par:/
1074	if (p_iFrMin==iKTABS2)
1075	{
1076	INDX=naf_maxx(iKTABS2m1, RTAB);
1077	IFR=((INDX+1)<=(iKTABS2/2))?INDX:INDX-iKTABS2;
1078	}
1079	else if (p_iFrMin>=0)
1080	{/* p_iFrMin et p_iFrMax positif */
1081	IFR=INDX=p_iFrMin+naf_maxx(p_iFrMax-p_iFrMin, RTAB+p_iFrMin);
1082	}
1083	else if (p_iFrMax<0)
1084	{/* p_iFrMin et p_iFrMax negatif */
1085	INDX=naf_maxx(p_iFrMax-p_iFrMin, RTAB+iKTABS2+p_iFrMin);
1086	IFR=INDX+p_iFrMin;
1087	}
1088	else
1089	{/* p_iFrMin et p_iFrMax de signe differents */
1090	int INDXNeg;
1091	/* frequence negative */
1092	INDXNeg=p_iFrMin+naf_maxx(-1-p_iFrMin, RTAB+iKTABS2+p_iFrMin);
1093	/* frequence positive */
1094	INDX=naf_maxx(p_iFrMax, RTAB);
1095	if (RTAB[INDX]<RTAB[INDXNeg+iKTABS2])
1096	{
1097	IFR=INDXNeg;
1098	}
1099	else
1100	{
1101	IFR=INDX;
1102	}
1103	}
1104	FR = IFRFREFO2IPAS;
1105	if (g_NAFVariable.IPRT==1)
1106	{
1107	fprintf(g_NAFVariable.NFPRT,"IFRMIN=%d IFRMAX=%d IFR=%d FR=%g RTAB=%g INDX=%d KTABS2=%d\n",p_iFrMin, p_iFrMax,
1108	IFR,FR, RTAB[INDX],INDX,iKTABS2);
1109	}
1110	SYSFREE(RTAB);
1111	return FR;
1112	/* v0.96 M. GASTINEAU 14/01/99 : fin de modification */
1113	}/* END SUBROUTINE FFTMAX*/
1114	#endif /* NAF_USE_OPTIMIZE */
1115	/*!
1116	!
1117	SUBROUTINE FOUR1(DATA,NN,ISIGN)*/
1118	/v0.96 M. GASTINEAU 14/12/98 : modification du calcul de la puissance /
1119	void naf_four1(double *DATA, int NN, int ISIGN)
1120	{
1121	/*
1122	!************** FOURIER RAPIDE DE NUMERICAL RECIPES
1123	implicit none
1124	! (DATA,NN,ISIGN)
1125	integer :: NN,ISIGN
1126	REAL (8) :: DATA(2*NN)
1127	!
1128	integer :: N,I,J,M,MMAX,ISTEP
1129	REAL (8) :: THETA,WPR,WPI,WR,WI,TEMPR,TEMPI,WTEMP
1130	*/
1131	int N,I,J,M,MMAX,ISTEP;
1132	double THETA,WPR,WPI,WR,WI,TEMPR,TEMPI,WTEMP;
1133	N=2*NN ;
1134	J=1;
1135	for(I=1;I<=N; I+=2)
1136	{
1137	/! en fait N est pair donc I<= N-1 /
1138	if(J>I)
1139	{
1140	TEMPR=DATA[J];
1141	TEMPI=DATA[J+1];
1142	DATA[J]=DATA[I];
1143	DATA[J+1]=DATA[I+1];
1144	DATA[I]=TEMPR;
1145	DATA[I+1]=TEMPI;
1146	}
1147	M=N/2;
1148	while ((M>=2) && (J>M))
1149	{
1150	J-=M; /J=J-M;/
1151	M >>=1; /M=M/2;/
1152	};
1153	J +=M; /J=J+M/
1154	}
1155	MMAX=2;
1156	while (N>MMAX)
1157	{
1158	ISTEP=2*MMAX;
1159	THETA=6.28318530717959E0/(ISIGN*MMAX);
1160	/v0.96 M. GASTINEAU 14/12/98 : modification /
1161	/WPR=-2.E0pow(sin(0.5E0THETA),2);//remplacee par:/
1162	WTEMP=sin(0.5E0*THETA);
1163	WPR=-2.E0WTEMPWTEMP;
1164	/v0.96 M. GASTINEAU 14/12/98 : fin modification /
1165	WPI=sin(THETA);
1166	WR=1.E0;
1167	WI=0.E0;
1168	for ( M=1; M<=MMAX; M+=2)
1169	{
1170	for (I=M; I<=N; I+=ISTEP)
1171	{
1172	J=I+MMAX;
1173	TEMPR=WRDATA[J]-WIDATA[J+1];
1174	TEMPI=WRDATA[J+1]+WIDATA[J];
1175	/!* TEMPR=SNGL[WR]DATA[J]-SNGL[WI]DATA[J+1]
1176	!** TEMPI=SNGL[WR]DATA[J+1]+SNGL[WI]DATA[J]*/
1177	DATA[J]=DATA[I]-TEMPR;
1178	DATA[J+1]=DATA[I+1]-TEMPI;
1179	DATA[I]=DATA[I]+TEMPR;
1180	DATA[I+1]=DATA[I+1]+TEMPI;
1181	}
1182	WTEMP=WR;
1183	WR=WRWPR-WIWPI+WR;
1184	WI=WIWPR+WTEMPWPI+WI;
1185	}
1186	MMAX=ISTEP;
1187	}
1188	}/* END SUBROUTINE FOUR1*/
1189
1190
1191	/* SUBROUTINE PUISS2(NT,N2)*/
1192	void naf_puiss2(int NT, int *N2)
1193	{
1194	/!******************************************
1195	! CALCULE LA PLUS GRANDE PUISSANCE DE 2
1196	! CONTENUE DANS NT
1197	!*******************************************
1198	implicit none
1199	! (NT,N2)
1200	INTEGER :: N2,NT
1201	!
1202	integer :: n
1203	! */
1204	int N;
1205	N=NT;
1206	if (N==0)
1207	{
1208	*N2=0;
1209	return;
1210	}
1211	*N2=1;
1212	while (N>=2)
1213	{
1214	N2 <<=1;/N2=N22/
1215	N >>=1; /N=N/2/
1216
1217	}
1218	}/* END SUBROUTINE PUISS2 */
1219
1220	/* SUBROUTINE MAXX(N,T,INDX)*/
1221	/!******************************************
1222	! CALCULE LE MAX ET L'INDICE DU MAX
1223	!*******************************************
1224	implicit none
1225	! (N,T,INDX)
1226	integer :: N,INDX
1227	REAL (8) :: T(0:N)
1228	!
1229	integer :: J
1230	REAL (8) VMAX
1231	!*/
1232	/* v0.96 M. GASTINEAU 12/01/99 : optimisation passage par retour et non par la pile */
1233	#if NAF_USE_OPTIMIZE==0
1234	void naf_maxx(int N, double T, int INDX)
1235	{
1236	int J;
1237	double VMAX;
1238	VMAX=T[0];
1239	*INDX=0;
1240	for(J=1; J<=N; J++)
1241	{
1242	if (T[J]>VMAX)
1243	{
1244	VMAX=T[J];
1245	*INDX=J;
1246	}
1247	}
1248	}/* END SUBROUTINE MAXX*/
1249	#else /remplacee par:/
1250	int naf_maxx(int N, double *T)
1251	{
1252	int J;
1253	double VMAX;
1254	int INDX=0;
1255	VMAX=T[0];
1256	for(J=1; J<=N; J++)
1257	{
1258	if (T[J]>VMAX)
1259	{
1260	VMAX=T[J];
1261	INDX=J;
1262	}
1263	}
1264	return INDX;
1265	}
1266	#endif /NAF_USE_OPTIMIZE/
1267	/* v0.96 M. GASTINEAU 12/01/99 : fin optimisation */
1268	/* END SUBROUTINE MAXX*/
1269
1270	/* SUBROUTINE MODTAB(N,T)*/
1271	void naf_modtab(int N, double *T)
1272	{
1273	/!***************************************
1274	! SUPPRESSION DE CERTAINS TERMES
1275	! A MODIFIER ULTERIEUREMENT (13/12/90)
1276	!****************************************
1277	implicit none
1278	! (N,T)
1279	integer :: N
1280	REAL (8) :: T(0:N)
1281	!
1282	integer NOK,I
1283	! */
1284	int NOK, I;
1285	NOK = 40;
1286	for(I = 0; I<=NOK; I++)
1287	{
1288	T[I] = 0;
1289	T[N-I] = 0;
1290	}
1291	} /* END SUBROUTINE MODTAB*/
1292
1293
1294	/* SUBROUTINE MODFRE(NUMFR,A,B)*/
1295	/v0.96 M. GASTINEAU 12/01/99 : modification dans le cas ou ICPLX=0 /
1296	/v0.97 M. GASTINEAU 27/05/99 : modification dans le cas ou ICPLX=0 et FS=0 /
1297	void naf_modfre(int NUMFR, double A, double B)
1298	{
1299	/*!-----------------------------------------------------------------------
1300	! PERMET DE MODIFIER UNE AMPLITUDE DEJA CALCULEE QUAND
1301	! ON RETROUVE LE MEME TERME, A TOL PRES DANS LE DEVELOPPEMENT
1302	!
1303	! NUMFR INDICE DE LA FREQUENCE EXISTANT DEJA
1304	! A,B PARTIES REELLES ET IMAGINAIRES DE L'AMPLITUDE DE LA MODIF
1305	! A APPORTER A L'AMPLITUDE DE LA FREQUENCE NUMFR
1306	!
1307	! AUTRES PARAMETRES TRANSMIS PAR LE COMMON/CGRAM/g_NAFVariable.ZAMP,g_NAFVariable.ZALP,g_NAFVariable.TFS,g_NAFVariable.NFS
1308	!
1309	! 30/4/91
1310	!-----------------------------------------------------------------------
1311
1312	!----------! g_NAFVariable.TFS EST LE TABLEAU FINAL DES FREQUENCES
1313	!----------! g_NAFVariable.ZAMP TABLEAU DES AMPLITUDES COMPLEXES
1314	!----------! g_NAFVariable.ZALP MATRICE DE PASSAGE DE LA BASE ORTHONORMALISEE
1315	!----------! g_NAFVariable.NFS NOMBRE DE FREQ.DEJA DETERMINEES
1316	IMPLICIT NONE
1317	! (NUMFR,A,B)
1318	integer :: NUMFR
1319	REAL (8) :: A,B
1320	!
1321	integer :: IT
1322	complex (8) :: ZI,ZOM,ZA,ZEX,ZINC
1323	complex (8), dimension (:), allocatable :: ZT
1324	! */
1325	int IT;
1326	t_complexe ZI,ZOM,ZA,ZEX,ZINC;
1327	t_complexe *ZT=NULL;
1328
1329	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
1330	#if NAF_USE_OPTIMIZE==0
1331	SYSCHECKMALLOCSIZE(ZT,t_complexe, g_NAFVariable.KTABS+1); /allocate(ZT(0:g_NAFVariable.KTABS),stat = NERROR)/
1332	ZI = cmplx(0.E0,1.E0);
1333	ZOM=muldoublcomplexe(g_NAFVariable.TFS[NUMFR]/g_NAFVariable.UNIANG,ZI); /ZOM=g_NAFVariable.TFS[NUMFR]/g_NAFVariable.UNIANGZI*/
1334	ZA=cmplx(A,B);
1335	fprintf(g_NAFVariable.NFPRT,"CORRECTION DE IFR = %d AMPLITUDE = %g\n",NUMFR, module(ZA));
1336	/*!-----------! L' AMPLITUDES DU TERMES EST CORRIGEES
1337	!-----------! ATTENTION ICI (CAS REEL) ON AURA AUSSI LE TERME CONJUGUE
1338	!-----------! QU'ON NE CALCULE PAS. LE TERME TOTAL EST
1339	!-----------! 2RE(g_NAFVariable.ZAMP(I)EXP(ZIg_NAFVariable.TFS(I)T) )*/
1340	g_NAFVariable.ZAMP[NUMFR]=addcomplexe(g_NAFVariable.ZAMP[NUMFR],ZA);
1341	if (g_NAFVariable.IPRT==1)
1342	{
1343	fprintf(g_NAFVariable.NFPRT," %+-20.15E %+-20.15E %+-20.15E %+-20.15E %+-20.15E\n",g_NAFVariable.TFS[NUMFR],
1344	module(g_NAFVariable.ZAMP[NUMFR]),g_NAFVariable.ZAMP[NUMFR].reel,
1345	g_NAFVariable.ZAMP[NUMFR].imag,atan2(g_NAFVariable.ZAMP[NUMFR].imag,g_NAFVariable.ZAMP[NUMFR].reel));
1346	}
1347	/!-----------! ON RETIRE LA CONTRIBUTION DU TERME g_NAFVariable.TFS(NUMFR) DANS TABS/
1348	if (g_NAFVariable.ICPLX==1)
1349	{
1350	ZEX=mulcomplexe(ZA,expcomplexe(muldoublcomplexe(g_NAFVariable.T0-g_NAFVariable.XH,ZOM))); /ZEX = ZAEXP(ZOM(g_NAFVariable.T0-g_NAFVariable.XH))/
1351	ZINC=expcomplexe(muldoublcomplexe(g_NAFVariable.XH,ZOM)); /ZINC=EXP(ZOMg_NAFVariable.XH)*/
1352	naf_ztpow(g_NAFVariable.KTABS,64,ZT,ZINC,ZEX); /CALL ZTPOW (KTABS+1,64,ZT,ZINC,ZEX)/
1353	for(IT=0;IT<=g_NAFVariable.KTABS;IT++)
1354	{
1355	g_NAFVariable.ZTABS[IT]=subcomplexe(g_NAFVariable.ZTABS[IT], ZT[IT]);
1356	}
1357	}
1358	else
1359	{
1360	ZEX=mulcomplexe(ZA,expcomplexe(muldoublcomplexe(g_NAFVariable.T0-g_NAFVariable.XH,ZOM))); /ZEX = ZAEXP(ZOM(g_NAFVariable.T0-g_NAFVariable.XH))/
1361	ZINC=expcomplexe(muldoublcomplexe(g_NAFVariable.XH,ZOM)); /ZINC=EXP(ZOMXH)*/
1362	naf_ztpow(g_NAFVariable.KTABS,64,ZT,ZINC,ZEX); /CALL ZTPOW (KTABS+1,64,ZT,ZINC,ZEX)/
1363	/v0.97 M. GASTINEAU 27/05/99 : ajout - correction bug si FS==0 /
1364	if (g_NAFVariable.TFS[NUMFR]==0.E0)
1365	{
1366	for(IT=0;IT<=g_NAFVariable.KTABS;IT++)
1367	{
1368	g_NAFVariable.ZTABS[IT].reel -= ZT[IT].reel;
1369	/g_NAFVariable.ZTABS(IT+1)=g_NAFVariable.ZTABS(IT+1)- DREAL(ZT(IT+1)) /
1370	}
1371	}
1372	else
1373	/v0.97 M. GASTINEAU 27/05/99 : fin ajout /
1374	for(IT=0;IT<=g_NAFVariable.KTABS;IT++)
1375	{
1376	/v0.96 M. GASTINEAU 12/01/99 : modification /
1377	/g_NAFVariable.ZTABS[IT].reel -=ZT[IT].reel;/
1378	/remplacee par:/
1379	g_NAFVariable.ZTABS[IT].reel -= 2*ZT[IT].reel;
1380	/v0.96 M. GASTINEAU 12/01/99 : fin modification /
1381	/* g_NAFVariable.ZTABS(IT)=g_NAFVariable.ZTABS(IT)- DREAL(ZT(IT)) */
1382	}
1383	}
1384	SYSFREE(ZT); /deallocate(ZT)/
1385	#else /remplacee par:/
1386	t_complexe pzarTabs, pzarZT;
1387	const int ikTabs=g_NAFVariable.KTABS; /v0.96 M. GASTINEAU 12/01/99 : optimisation/
1388
1389	SYSCHECKMALLOCSIZE(ZT,t_complexe, g_NAFVariable.KTABS+1); /allocate(ZT(0:g_NAFVariable.KTABS),stat = NERROR)/
1390	i_compl_cmplx(&ZI,0.E0,1.E0);
1391	ZOM=i_compl_muldoubl(g_NAFVariable.TFS[NUMFR]/g_NAFVariable.UNIANG,ZI); /ZOM=g_NAFVariable.TFS[NUMFR]/g_NAFVariable.UNIANGZI*/
1392
1393	i_compl_cmplx(&ZA,A,B);
1394	fprintf(g_NAFVariable.NFPRT,"CORRECTION DE IFR = %d AMPLITUDE = %g\n",NUMFR, i_compl_module(ZA));
1395	/*!-----------! L' AMPLITUDES DU TERMES EST CORRIGEES
1396	!-----------! ATTENTION ICI (CAS REEL) ON AURA AUSSI LE TERME CONJUGUE
1397	!-----------! QU'ON NE CALCULE PAS. LE TERME TOTAL EST
1398	!-----------! 2RE(g_NAFVariable.ZAMP(I)EXP(ZIg_NAFVariable.TFS(I)T) )*/
1399	i_compl_padd(g_NAFVariable.ZAMP+NUMFR,&ZA);
1400	if (g_NAFVariable.IPRT==1)
1401	{
1402	fprintf(g_NAFVariable.NFPRT," %+-20.15E %+-20.15E %+-20.15E %+-20.15E %+-20.15E\n",g_NAFVariable.TFS[NUMFR],
1403	i_compl_module(g_NAFVariable.ZAMP[NUMFR]),g_NAFVariable.ZAMP[NUMFR].reel,
1404	g_NAFVariable.ZAMP[NUMFR].imag,atan2(g_NAFVariable.ZAMP[NUMFR].imag,g_NAFVariable.ZAMP[NUMFR].reel));
1405	}
1406	/!-----------! ON RETIRE LA CONTRIBUTION DU TERME g_NAFVariable.TFS(NUMFR) DANS TABS/
1407	if (g_NAFVariable.ICPLX==1)
1408	{
1409	ZEX=i_compl_mul(ZA,i_compl_exp(i_compl_muldoubl(g_NAFVariable.T0-g_NAFVariable.XH,ZOM))); /ZEX = ZAEXP(ZOM(g_NAFVariable.T0-g_NAFVariable.XH))/
1410	ZINC=i_compl_exp(i_compl_muldoubl(g_NAFVariable.XH,ZOM)); /ZINC=EXP(ZOMg_NAFVariable.XH)*/
1411	naf_ztpow(g_NAFVariable.KTABS,64,ZT,ZINC,ZEX); /CALL ZTPOW (KTABS+1,64,ZT,ZINC,ZEX)/
1412	/v0.96 M. GASTINEAU 12/01/99 : optimisation/
1413	/*for(IT=0, pzarTabs=g_NAFVariable.ZTABS, pzarZT = ZT;
1414	IT<=g_NAFVariable.KTABS;
1415	IT++, pzarTabs++, pzarZT++)//remplacee par:*/
1416	for(IT=0, pzarTabs=g_NAFVariable.ZTABS, pzarZT = ZT;
1417	IT<=ikTabs;
1418	IT++, pzarTabs++, pzarZT++)/v0.96 M. GASTINEAU 12/01/99 : fin optimisation/
1419	{
1420	i_compl_psub(pzarTabs,pzarZT);
1421	}
1422	}
1423	else
1424	{
1425	ZEX=i_compl_mul(ZA,i_compl_exp(i_compl_muldoubl(g_NAFVariable.T0-g_NAFVariable.XH,ZOM))); /ZEX = ZAEXP(ZOM(g_NAFVariable.T0-g_NAFVariable.XH))/
1426	ZINC=i_compl_exp(i_compl_muldoubl(g_NAFVariable.XH,ZOM)); /ZINC=EXP(ZOMXH)*/
1427	naf_ztpow(g_NAFVariable.KTABS,64,ZT,ZINC,ZEX); /CALL ZTPOW (KTABS+1,64,ZT,ZINC,ZEX)/
1428	/v0.96 M. GASTINEAU 12/01/99 : optimisation/
1429	/for(IT=0;IT<=g_NAFVariable.KTABS;IT++)/
1430	/v0.97 M. GASTINEAU 27/05/99 : ajout - correction bug si FS==0 /
1431	if (g_NAFVariable.TFS[NUMFR]==0.E0)
1432	{
1433	for(IT=0;IT<=ikTabs;IT++)
1434	{
1435	g_NAFVariable.ZTABS[IT].reel -= ZT[IT].reel;
1436	/g_NAFVariable.ZTABS(IT+1)=g_NAFVariable.ZTABS(IT+1)- DREAL(ZT(IT+1)) /
1437	}
1438	}
1439	else
1440	/v0.97 M. GASTINEAU 27/05/99 : fin ajout /
1441	for(IT=0;IT<=ikTabs;IT++)/v0.96 M. GASTINEAU 12/01/99 : fin optimisation/
1442	{
1443	/v0.96 M. GASTINEAU 12/01/99 : modification /
1444	/* g_NAFVariable.ZTABS[IT].reel -=ZT[IT].reel;// g_NAFVariable.ZTABS(IT)=g_NAFVariable.ZTABS(IT)- DREAL(ZT(IT)) */
1445	/remplacee par:/
1446	g_NAFVariable.ZTABS[IT].reel -=2ZT[IT].reel;/ g_NAFVariable.ZTABS(IT)=g_NAFVariable.ZTABS(IT)- DREAL(ZT(IT)) */
1447	/v0.96 M. GASTINEAU 12/01/99 : fin modification /
1448	}
1449	}
1450	SYSFREE(ZT); /deallocate(ZT)/
1451	#endif /NAF_USE_OPTIMIZE==0/
1452	/* v0.96 M. GASTINEAU 01/12/98 : fin optimisation */
1453	}/* END SUBROUTINE MODFRE*/
1454
1455
1456	/* SUBROUTINE GRAMSC(FS,A,B)*/
1457	/v0.96 M. GASTINEAU 12/01/99 : modification dans le cas ou ICPLX=0 /
1458	/v0.97 M. GASTINEAU 27/05/99 : modification dans le cas ou ICPLX=0 et FS=0 /
1459	BOOL naf_gramsc(double FS, double A, double B)
1460	{
1461	/*
1462	!-----------------------------------------------------------------------
1463	! GRAMSC ORTHONORMALISATION PAR GRAM-SCHIMDT DE LA BASE DE FONCTI
1464	! CONSTITUEE DES TERMES PERIODIQUES TROUVES PAR ANALYSE
1465	! SPECTRALE .
1466	! FS FREQUENCE EN "/AN
1467	! A,B PARTIES REELLES ET IMAGINAIRES DE L'AMPLITUDE
1468	!
1469	! AUTRES PARAMETRES TRANSMIS PAR LE COMMON/CGRAM/g_NAFVariable.ZAMP,g_NAFVariable.ZALP,g_NAFVariable.TFS,g_NAFVariable.NFS
1470	!
1471	! MODIFIE LE 26 9 87 POUR LES FONCTIONS REELLES J. LASKAR
1472	!-----------------------------------------------------------------------
1473
1474	!----------! g_NAFVariable.TFS EST LE TABLEAU FINAL DES FREQUENCES
1475	!----------! g_NAFVariable.ZAMP TABLEAU DES AMPLITUDES COMPLEXES
1476	!----------! g_NAFVariable.ZALP MATRICE DE PASSAGE DE LA BASE ORTHONORMALISEE
1477	!----------! ZTEE TABLEAU DE TRAVAIL
1478	!----------! g_NAFVariable.NFS NOMBRE DE FREQ.DEJA DETERMINEES
1479	IMPLICIT NONE
1480	! (g_NAFVariable.KTABS,FS,A,B)
1481	REAL (8) :: FS,A,B
1482	!
1483	integer :: I, J, K, NF, IT
1484	REAL (8) DIV
1485	complex (8), dimension(:),allocatable :: ZTEE
1486	complex (8) :: ZDIV,ZMUL,ZI,ZEX,ZINC,ZA,ZOM
1487	complex (8), dimension (:),allocatable :: ZT
1488	! */
1489	int I, J, K, NF, IT;
1490	double DIV;
1491	t_complexe *ZTEE=NULL;
1492	t_complexe ZDIV,ZMUL,ZI,ZEX,ZINC,ZA,ZOM;
1493	t_complexe *ZT=NULL;
1494	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
1495	#if NAF_USE_OPTIMIZE>0
1496	t_complexe pzarTabs, pzarZT;
1497	const int ikTabs=g_NAFVariable.KTABS; /v0.96 M. GASTINEAU 12/01/99 : optimisation/
1498	#endif /**/
1499	SYSCHECKMALLOCSIZE(ZT, t_complexe, g_NAFVariable.KTABS+1); /allocate(ZT(0:g_NAFVariable.KTABS),stat = NERROR)/
1500	SYSCHECKMALLOCSIZE(ZTEE, t_complexe, g_NAFVariable.NTERM+1); /allocate(ZTEE(1:g_NAFVariable.NTERM),stat = NERROR)/
1501
1502	/!----------! CALCUL DE ZTEE(I)=<EN,EI>/
1503	for(I =1;I<=g_NAFVariable.NFS;I++)
1504	{
1505	if(naf_proscaa(FS,g_NAFVariable.TFS[I],ZTEE+I)==FALSE)
1506	{
1507	SYSFREE(ZT);
1508	SYSFREE(ZTEE);
1509	return FALSE;
1510	}
1511	}
1512	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
1513	#if NAF_USE_OPTIMIZE==0
1514	/!----------! NF EST LE NUMERO DE LA NOUVELLE FREQUENCE/
1515	NF=g_NAFVariable.NFS+1;
1516	ZTEE[NF]=cmplx(1.E0,0.E0);
1517	/!----------! CALCUL DE FN = EN - SOM(<EN,FI>FI) QUI EST ORTHOGONAL AUX F/
1518	g_NAFVariable.TFS[NF]=FS;
1519	for( K=1;K<=g_NAFVariable.NFS;K++)
1520	{
1521	for(I=K;I<=g_NAFVariable.NFS;I++)
1522	{
1523	for(J=1;J<=I;J++)
1524	{
1525	/g_NAFVariable.ZALP(NF,K)=g_NAFVariable.ZALP(NF,K)-DCONJG(g_NAFVariable.ZALP(I,J))g_NAFVariable.ZALP(I,K)ZTEE(J)/
1526	g_NAFVariable.ZALP[NF][K]=subcomplexe(g_NAFVariable.ZALP[NF][K],mulcomplexe(mulcomplexe(conjcomplexe(g_NAFVariable.ZALP[I][J]),g_NAFVariable.ZALP[I][K]),ZTEE[J]));
1527	}
1528	}
1529	}
1530	g_NAFVariable.ZALP[NF][NF]=cmplx(1.E0,0.E0);
1531	/!----------! ON REND LA NORME DE FN = 1/
1532	DIV=1.E0;
1533	ZDIV=cmplx(0.E0,0.E0);
1534	for( I=1; I<=NF; I++)
1535	{
1536	/ZDIV=ZDIV+DCONJG(g_NAFVariable.ZALP(NF,I))ZTEE(I)*/
1537	ZDIV=addcomplexe(ZDIV, mulcomplexe(conjcomplexe(g_NAFVariable.ZALP[NF][I]),ZTEE[I]));
1538	}
1539	DIV=sqrt(module(ZDIV));
1540	if (g_NAFVariable.IPRT==1)
1541	{
1542	/v0.96 M. GASTINEAU 19/11/98 : modification/
1543	/fprintf(g_NAFVariable.NFPRT,"ZDIV= %g+i%g DIV=%g\n",ZDIV,DIV);//remplacee par:/
1544	fprintf(g_NAFVariable.NFPRT,"ZDIV= %g+i%g DIV=%g\n",ZDIV.reel,ZDIV.imag,DIV);
1545	/v0.96 M. GASTINEAU 19/11/98 : fin modification/
1546	}
1547	for(I=1; I<=NF; I++)
1548	{
1549	g_NAFVariable.ZALP[NF][I]=muldoublcomplexe(1.E0/DIV,g_NAFVariable.ZALP[NF][I]); /g_NAFVariable.ZALP(NF,I) = g_NAFVariable.ZALP(NF,I)/DIV/
1550	}
1551	/*!-----------! F1,F2....., FN EST UNE BASE ORTHONORMEE
1552	!-----------! ON RETIRE MAINTENANT A F <F,FN>FN (<F,FN>=<F,EN>)*/
1553	ZMUL=muldoublcomplexe(1.E0/DIV, cmplx(A,B));
1554	ZI=cmplx(0.E0,1.E0);
1555	for(I=1; I<=NF; I++)
1556	{
1557	ZOM=muldoublcomplexe(g_NAFVariable.TFS[I]/g_NAFVariable.UNIANG,ZI); /ZOM=g_NAFVariable.TFS(I)/g_NAFVariable.UNIANGZI*/
1558	ZA=mulcomplexe(g_NAFVariable.ZALP[NF][I],ZMUL);
1559	/*!-----------! LES AMPLITUDES DES TERMES SONT CORRIGEES
1560	!-----------! ATTENTION ICI (CAS REEL) ON AURA AUSSI LE TERME CONJUGUE
1561	!-----------! QU'ON NE CALCULE PAS. LE TERME TOTAL EST
1562	!-----------! 2RE(g_NAFVariable.ZAMP(I)EXP(ZIg_NAFVariable.TFS(I)T) )*/
1563	g_NAFVariable.ZAMP[I]=addcomplexe(g_NAFVariable.ZAMP[I],ZA);
1564	if (g_NAFVariable.IPRT==1)
1565	{
1566	fprintf(g_NAFVariable.NFPRT," %g %g %g %g %g\n", g_NAFVariable.TFS[I],module(g_NAFVariable.ZAMP[I]),g_NAFVariable.ZAMP[I].reel,
1567	g_NAFVariable.ZAMP[I].imag,atan2(g_NAFVariable.ZAMP[I].imag,g_NAFVariable.ZAMP[I].reel));
1568	}
1569	/!-----------! ON RETIRE LA CONTRIBUTION TOTALE DU TERME g_NAFVariable.TFS(I) DANS TABS/
1570	if (g_NAFVariable.ICPLX==1)
1571	{
1572	ZEX=mulcomplexe(ZA, expcomplexe(muldoublcomplexe(g_NAFVariable.T0-g_NAFVariable.XH,ZOM))); /ZEX = ZAEXP(ZOM(g_NAFVariable.T0-g_NAFVariable.XH))/
1573	ZINC=expcomplexe(muldoublcomplexe(g_NAFVariable.XH,ZOM)); /ZINC=EXP(ZOMg_NAFVariable.XH)*/
1574	naf_ztpow(g_NAFVariable.KTABS,64,ZT,ZINC,ZEX); /naf_ztpow(g_NAFVariable.KTABS+1,64,ZT,ZINC,ZEX);/
1575	for(IT=0;IT<=g_NAFVariable.KTABS;IT++)
1576	{
1577	g_NAFVariable.ZTABS[IT]=subcomplexe(g_NAFVariable.ZTABS[IT],ZT[IT]);
1578	}
1579	}
1580	else
1581	{
1582	ZEX=mulcomplexe(ZA, expcomplexe(muldoublcomplexe(g_NAFVariable.T0-g_NAFVariable.XH,ZOM))); /ZEX = ZAEXP(ZOM(g_NAFVariable.T0-g_NAFVariable.XH))/
1583	ZINC=expcomplexe(muldoublcomplexe(g_NAFVariable.XH,ZOM)); /ZINC=EXP(ZOMg_NAFVariable.XH)*/
1584	naf_ztpow(g_NAFVariable.KTABS,64,ZT,ZINC,ZEX);/naf_ztpow(g_NAFVariable.KTABS+1,64,ZT,ZINC,ZEX);/
1585	/v0.97 M. GASTINEAU 27/05/99 : ajout - correction bug si FS==0 /
1586	if (FS==0.E0)
1587	{
1588	for(IT=0;IT<=g_NAFVariable.KTABS;IT++)
1589	{
1590	g_NAFVariable.ZTABS[IT].reel -= ZT[IT].reel;
1591	/g_NAFVariable.ZTABS(IT+1)=g_NAFVariable.ZTABS(IT+1)- DREAL(ZT(IT+1)) /
1592	}
1593	}
1594	else
1595	/v0.97 M. GASTINEAU 27/05/99 : fin ajout /
1596	for(IT=0;IT<=g_NAFVariable.KTABS;IT++)
1597	{
1598	/v0.96 M. GASTINEAU 12/01/99 : modification/
1599	/g_NAFVariable.ZTABS[IT].reel -=ZT[IT].reel;/ /g_NAFVariable.ZTABS(IT+1)=g_NAFVariable.ZTABS(IT+1)- DREAL(ZT(IT+1)) /
1600	/remplacee par: /
1601	g_NAFVariable.ZTABS[IT].reel -= 2ZT[IT].reel; /g_NAFVariable.ZTABS(IT+1)=g_NAFVariable.ZTABS(IT+1)- DREAL(ZT(IT+1)) */
1602	/v0.96 M. GASTINEAU 12/01/99 : fin modification/
1603	}
1604	}
1605	}
1606	#else /remplacee par:/
1607	/!----------! NF EST LE NUMERO DE LA NOUVELLE FREQUENCE/
1608	NF=g_NAFVariable.NFS+1;
1609	i_compl_cmplx(ZTEE+NF,1.E0,0.E0);
1610	/!----------! CALCUL DE FN = EN - SOM(<EN,FI>FI) QUI EST ORTHOGONAL AUX F/
1611	g_NAFVariable.TFS[NF]=FS;
1612	for( K=1;K<=g_NAFVariable.NFS;K++)
1613	{
1614	for(I=K;I<=g_NAFVariable.NFS;I++)
1615	{
1616	for(J=1;J<=I;J++)
1617	{
1618	/g_NAFVariable.ZALP(NF,K)=g_NAFVariable.ZALP(NF,K)-DCONJG(g_NAFVariable.ZALP(I,J))g_NAFVariable.ZALP(I,K)ZTEE(J)/
1619	t_complexe zSubExp;
1620	zSubExp=i_compl_mul(i_compl_mul(i_compl_conj(&(g_NAFVariable.ZALP[I][J])),g_NAFVariable.ZALP[I][K]),ZTEE[J]);
1621	i_compl_psub(&(g_NAFVariable.ZALP[NF][K]),&zSubExp);
1622	}
1623	}
1624	}
1625	i_compl_cmplx(&(g_NAFVariable.ZALP[NF][NF]),1.E0,0.E0);
1626	/!----------! ON REND LA NORME DE FN = 1/
1627	DIV=1.E0;
1628	i_compl_cmplx(&ZDIV,0.E0,0.E0);
1629	for( I=1; I<=NF; I++)
1630	{
1631	/ZDIV=ZDIV+DCONJG(g_NAFVariable.ZALP(NF,I))ZTEE(I)*/
1632	t_complexe zSubExp;
1633	zSubExp = i_compl_mul(i_compl_conj(&(g_NAFVariable.ZALP[NF][I])),ZTEE[I]);
1634	i_compl_padd(&ZDIV,&zSubExp);
1635	}
1636	DIV=sqrt(i_compl_module(ZDIV));
1637	if (g_NAFVariable.IPRT==1)
1638	{
1639	/v0.96 M. GASTINEAU 19/11/98 : modification/
1640	/fprintf(g_NAFVariable.NFPRT,"ZDIV= %g+i%g DIV=%g\n",ZDIV,DIV);//remplacee par:/
1641	fprintf(g_NAFVariable.NFPRT,"ZDIV= %g+i%g DIV=%g\n",ZDIV.reel,ZDIV.imag,DIV);
1642	/v0.96 M. GASTINEAU 19/11/98 : fin modification/
1643	}
1644	for(I=1; I<=NF; I++)
1645	{
1646	i_compl_pdivdoubl(&(g_NAFVariable.ZALP[NF][I]),&DIV); /g_NAFVariable.ZALP(NF,I) = g_NAFVariable.ZALP(NF,I)/DIV/
1647	}
1648	/*!-----------! F1,F2....., FN EST UNE BASE ORTHONORMEE
1649	!-----------! ON RETIRE MAINTENANT A F <F,FN>FN (<F,FN>=<F,EN>)*/
1650	/v0.96 M. GASTINEAU 01/12/98 : modification /
1651	i_compl_cmplx(&ZMUL,A/DIV,B/DIV);
1652	i_compl_cmplx(&ZI,0.E0,1.E0);
1653	for(I=1; I<=NF; I++)
1654	{
1655	ZOM=i_compl_muldoubl(g_NAFVariable.TFS[I]/g_NAFVariable.UNIANG,ZI); /ZOM=g_NAFVariable.TFS(I)/g_NAFVariable.UNIANGZI*/
1656	ZA=i_compl_mul(g_NAFVariable.ZALP[NF][I],ZMUL);
1657	/*!-----------! LES AMPLITUDES DES TERMES SONT CORRIGEES
1658	!-----------! ATTENTION ICI (CAS REEL) ON AURA AUSSI LE TERME CONJUGUE
1659	!-----------! QU'ON NE CALCULE PAS. LE TERME TOTAL EST
1660	!-----------! 2RE(g_NAFVariable.ZAMP(I)EXP(ZIg_NAFVariable.TFS(I)T) )*/
1661	i_compl_padd(g_NAFVariable.ZAMP+I,&ZA);
1662	if (g_NAFVariable.IPRT==1)
1663	{
1664	fprintf(g_NAFVariable.NFPRT," %g %g %g %g %g\n", g_NAFVariable.TFS[I],i_compl_module(g_NAFVariable.ZAMP[I]),g_NAFVariable.ZAMP[I].reel,
1665	g_NAFVariable.ZAMP[I].imag,atan2(g_NAFVariable.ZAMP[I].imag,g_NAFVariable.ZAMP[I].reel));
1666	}
1667	/!-----------! ON RETIRE LA CONTRIBUTION TOTALE DU TERME g_NAFVariable.TFS(I) DANS TABS/
1668	if (g_NAFVariable.ICPLX==1)
1669	{
1670	ZEX=i_compl_mul(ZA, i_compl_exp(i_compl_muldoubl(g_NAFVariable.T0-g_NAFVariable.XH,ZOM))); /ZEX = ZAEXP(ZOM(g_NAFVariable.T0-g_NAFVariable.XH))/
1671	ZINC=i_compl_exp(i_compl_muldoubl(g_NAFVariable.XH,ZOM)); /ZINC=EXP(ZOMg_NAFVariable.XH)*/
1672	naf_ztpow(g_NAFVariable.KTABS,64,ZT,ZINC,ZEX); /naf_ztpow(g_NAFVariable.KTABS+1,64,ZT,ZINC,ZEX);/
1673	/v0.96 M. GASTINEAU 12/01/99 : optimisation/
1674	/*for(IT=0, pzarTabs=g_NAFVariable.ZTABS, pzarZT=ZT;
1675	IT<=g_NAFVariable.KTABS;
1676	IT++,pzarTabs++,pzarZT++)//remplacee par:*/
1677	for(IT=0, pzarTabs=g_NAFVariable.ZTABS, pzarZT=ZT;
1678	IT<=ikTabs;
1679	IT++,pzarTabs++,pzarZT++) /v0.96 M. GASTINEAU 12/01/99 : optimisation/
1680	{
1681	i_compl_psub(pzarTabs, pzarZT);
1682	}
1683	}
1684	else
1685	{
1686	ZEX=i_compl_mul(ZA, i_compl_exp(i_compl_muldoubl(g_NAFVariable.T0-g_NAFVariable.XH,ZOM))); /ZEX = ZAEXP(ZOM(g_NAFVariable.T0-g_NAFVariable.XH))/
1687	ZINC=i_compl_exp(i_compl_muldoubl(g_NAFVariable.XH,ZOM)); /ZINC=EXP(ZOMg_NAFVariable.XH)*/
1688	naf_ztpow(g_NAFVariable.KTABS,64,ZT,ZINC,ZEX);/naf_ztpow(g_NAFVariable.KTABS+1,64,ZT,ZINC,ZEX);/
1689	/v0.96 M. GASTINEAU 12/01/99 : optimisation/
1690	/for(IT=0;IT<=g_NAFVariable.KTABS;IT++)//remplacee par:/
1691	/v0.97 M. GASTINEAU 27/05/99 : ajout - correction bug si FS==0 /
1692	if (FS==0.E0)
1693	{
1694	for(IT=0;IT<=ikTabs;IT++)
1695	{
1696	g_NAFVariable.ZTABS[IT].reel -= ZT[IT].reel;
1697	/g_NAFVariable.ZTABS(IT+1)=g_NAFVariable.ZTABS(IT+1)- DREAL(ZT(IT+1)) /
1698	}
1699	}
1700	else
1701	/v0.97 M. GASTINEAU 27/05/99 : fin ajout /
1702	for(IT=0;IT<=ikTabs;IT++)
1703	/v0.96 M. GASTINEAU 12/01/99 : fin optimisation/
1704	{
1705	/v0.96 M. GASTINEAU 12/01/99 : modification/
1706	/* g_NAFVariable.ZTABS[IT].reel -=ZT[IT].reel; //g_NAFVariable.ZTABS(IT+1)=g_NAFVariable.ZTABS(IT+1)- DREAL(ZT(IT+1)) */
1707	/remplacee par:/
1708	g_NAFVariable.ZTABS[IT].reel -= 2ZT[IT].reel; /g_NAFVariable.ZTABS(IT+1)=g_NAFVariable.ZTABS(IT+1)- DREAL(ZT(IT+1)) */
1709	/v0.96 M. GASTINEAU 12/01/99 : fin modification/
1710	}
1711	}
1712	}
1713	#endif /NAF_USE_OPTIMIZE==0/
1714	/* v0.96 M. GASTINEAU 01/12/98 : fin optimisation */
1715	SYSFREE(ZT);
1716	SYSFREE(ZTEE);
1717	return TRUE;
1718	}/* END SUBROUTINE GRAMSC*/
1719
1720
1721	/* SUBROUTINE PROSCA(F1,F2,ZP)*/
1722	void naf_prosca(double F1, double F2, t_complexe* ZP)
1723	{
1724	/*!-----------------------------------------------------------------------
1725	! PROSCA CALCULE LE PRODUIT SCALAIRE DE EXP(IF1T) PAR EXP (-I*F
1726	! SUR L'INTERVALLE [0:g_NAFVariable.KTABS] T=g_NAFVariable.T0+g_NAFVariable.XH*IT
1727	! CALCUL ANALYTIQUE
1728	! ZP PRODUIT SCALAIRE COMPLEXE
1729	! F1,F2 FREQUENCES EN "/AN
1730	! REVU LE 26/9/87 J. LASKAR
1731	!-----------------------------------------------------------------------
1732
1733	IMPLICIT NONE
1734	! (g_NAFVariable.KTABS,F1,F2,ZP)
1735	REAL (8) :: F1,F2
1736	complex (8) :: ZP
1737	!
1738	complex (8) ZI,ZF,ZF1,ZF2
1739	REAL (8) PICARRE, T1,T2, XT,DIV,FAC,T,FR1,FR2
1740	* */
1741	t_complexe ZI,ZF,ZF1,ZF2;
1742	double PICARRE, T1,T2, XT,DIV,FACTEUR,T,FR1,FR2;
1743	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
1744	#if NAF_USE_OPTIMIZE==0
1745	ZI=cmplx(0.E0,1.E0);
1746	#else /remplacee par:/
1747	i_compl_cmplx(&ZI,0.E0,1.E0);
1748	#endif /NAF_USE_OPTIMIZE==0/
1749	/* v0.96 M. GASTINEAU 01/12/98 : fin modification*/
1750	/!----------! FREQUENCES EN UNITE D'ANGLE PAR UNITE DE TEMPS/
1751	FR1=F1/g_NAFVariable.UNIANG;
1752	FR2=F2/g_NAFVariable.UNIANG;
1753	T1 =g_NAFVariable.T0;
1754	T2 =g_NAFVariable.T0+g_NAFVariable.KTABS*g_NAFVariable.XH;
1755	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
1756	#if NAF_USE_OPTIMIZE==0
1757	/v0.96 M. GASTINEAU 01/12/98 : correction bug /
1758	/ZF=muldoublcomplexe((FR1-FR2)(T2-T1), ZF);/ /ZF=ZI(FR1-FR2)(T2-T1)*/
1759	/remplacee par: /ZF=muldoublcomplexe((FR1-FR2)*(T2-T1), ZI);
1760	ZF1=muldoublcomplexe((FR1-FR2)T1, ZI); /ZF1=ZI(FR1-FR2)T1*/
1761	ZF2=muldoublcomplexe((FR1-FR2)T2, ZI); /ZF2=ZI(FR1-FR2)T2*/
1762	ZP=divcomplexe(subcomplexe(expcomplexe(ZF2),expcomplexe(ZF1)),ZF); /ZP=(EXP(ZF2)-EXP(ZF1))/ZF*/
1763	#else /remplacee par:/
1764	ZF=i_compl_muldoubl((FR1-FR2)(T2-T1), ZI); /ZF=ZI(FR1-FR2)(T2-T1)*/
1765	ZF1=i_compl_muldoubl((FR1-FR2)T1, ZI); /ZF1=ZI(FR1-FR2)T1*/
1766	ZF2=i_compl_muldoubl((FR1-FR2)T2, ZI); /ZF2=ZI(FR1-FR2)T2*/
1767	ZP=i_compl_div(i_compl_sub(i_compl_exp(ZF2),i_compl_exp(ZF1)),ZF); /ZP=(EXP(ZF2)-EXP(ZF1))/ZF*/
1768	#endif /NAF_USE_OPTIMIZE==0/
1769	/* v0.96 M. GASTINEAU 01/12/98 : fin optimisation */
1770	/*!
1771	! PI=2.D0ATAN2(1.D0,0.D0)/
1772	PICARRE=pi*pi;
1773	T=(T2-T1)/2.E0;
1774	XT=(FR1-FR2)*T;
1775	DIV=XT*XT-PICARRE;
1776	if(fabs(DIV)<1.E-10)
1777	{
1778	/v0.96 M. GASTINEAU 19/11/98 : modification/
1779	/*fprintf(stdout, "g_NAFVariable.T0= %g, g_NAFVariable.XH=%g, g_NAFVariable.KTABS==%g\n",
1780	g_NAFVariable.T0,g_NAFVariable.XH,g_NAFVariable.KTABS);*/
1781	/remplacee par:/
1782	fprintf(stdout, "g_NAFVariable.T0= %g, g_NAFVariable.XH=%g, g_NAFVariable.KTABS=%d\n",
1783	g_NAFVariable.T0, g_NAFVariable.XH, g_NAFVariable.KTABS);
1784	/v0.96 M. GASTINEAU 19/11/98 : fin modification/
1785	fprintf(stdout, "F1=%g ,F2=%g\n",F1,F2);
1786	fprintf(stdout, "T1= %g , T2= %g, T= %g\n",T1, T2, T);
1787	fprintf(stdout, "FR1= %g ,FR2= %g ,XT= %g \n",FR1,FR2,XT);
1788	fprintf(stdout, "DIV= %g\n",DIV);
1789	} else
1790	{
1791	FACTEUR=-PICARRE/DIV;
1792	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
1793	#if NAF_USE_OPTIMIZE==0
1794	ZP=muldoublcomplexe(FACTEUR,ZP);
1795	#else /remplacee par:/
1796	i_compl_pmuldoubl(ZP,&FACTEUR);
1797	#endif /NAF_USE_OPTIMIZE==0/
1798	/* v0.96 M. GASTINEAU 01/12/98 : fin optimisation */
1799	}
1800	}/* END SUBROUTINE PROSCA*/
1801
1802	/* FUNCTION FUNC(X) */
1803	double naf_func(double X)
1804	{
1805	/* IMPLICIT NONE
1806	! (X)
1807	REAL (8) :: X,FUNC
1808	!
1809	REAL (8) ::RMD
1810	! */
1811	double RMD = 0.0;
1812	naf_profre(X,&AF,&BF,&RMD);
1813	return RMD;
1814	}/* END FUNCTION FUNC*/
1815
1816	/* FUNCTION FUNCP(X)*/
1817	double naf_funcp(double X)
1818	{
1819	/* IMPLICIT NONE
1820	! (X)
1821	REAL (8) :: X,FUNCP
1822	!
1823	REAL (8) :: DER, RMF
1824	! */
1825	double DER, RMF;
1826	naf_proder(X,&DER,&AF,&BF,&RMF);
1827	return DER;
1828	}/* END FUNCTION FUNCP*/
1829
1830	/* SUBROUTINE PRODER(FS,DER,A,B,RM)*/
1831	void naf_proder(double FS, double DER, double A, double B,double RM)
1832	{
1833	/!**********************************************************************
1834	! CE SOUS PROG. CALCULE LA DERIVEE DU CARRE DU MODULE DU PRODUIT
1835	! SCALAIRE DE EXP(-IFS)T PAR TAB(0:7
1836	! UTILISE LA METHODE D'INTEGRATION DE HARDY
1837	! LE RESULTAT EST DANS DER
1838	! A +IB LE PRODUIT SCALAIRE , RM SON MODULE
1839	! FS EST DONNEE EN SECONDES PAR AN
1840	! REVU LE 26/9/87 J. LASKAR
1841	! 26/5/98 (CHGT DE SIGNE) (*m/4)
1842	!***********************************************************************
1843
1844	IMPLICIT NONE
1845	! (FS,DER,A,B,RM)
1846	REAL (8) :: FS,DER,A,B,RM
1847	!
1848	integer,parameter :: NVECT=64
1849	REAL (8) :: OM,ANG0,ANGI,H
1850	complex (8) :: ZI,ZAC,ZINC,ZEX,ZB,ZT,ZA
1851	integer :: LTF
1852	complex (8), dimension (:), allocatable :: ZTF
1853	! */
1854	#define NVECT 64
1855	double OM,ANG0,ANGI,H;
1856	t_complexe ZI,ZAC,ZINC,ZEX,ZB,/ZT,/ZA;
1857	int LTF;
1858	t_complexe ZTF=NULL;/tableau de 1 a KTABS+1 */
1859	SYSCHECKMALLOCSIZE(ZTF,t_complexe, g_NAFVariable.KTABS+1); /allocate(ZTF(0:g_NAFVariable.KTABS),stat = NERROR)/
1860	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
1861	#if NAF_USE_OPTIMIZE==0
1862	/*!
1863	!------------------ CONVERSION DE FS EN RD/AN*/
1864	OM=FS/g_NAFVariable.UNIANG;
1865	ZI=cmplx(0.E0,1.E0);
1866	/*!------------------ LONGUEUR DU TABLEAU A INTEGRER
1867	!------------------ DANS LE CAS REEL MEME CHOSE */
1868	LTF=g_NAFVariable.KTABS;
1869	ANG0=OM*g_NAFVariable.T0;
1870	ANGI=OM*g_NAFVariable.XH;
1871	ZAC = expcomplexe (muldoublcomplexe(-ANG0,ZI));
1872	ZINC= expcomplexe (muldoublcomplexe(-ANGI,ZI));
1873	ZEX = divcomplexe(ZAC,ZINC); /ZEX = ZAC/ZINC/
1874	naf_ztpow2(g_NAFVariable.KTABS,NVECT,ZTF,g_NAFVariable.ZTABS,TWIN,ZINC,ZEX); /CALL ZTPOW2(g_NAFVariable.KTABS+1,NVECT,ZTF,g_NAFVariable.ZTABS,TWIN,ZINC,ZEX)/
1875	/!------------------ TAILLE DU PAS/
1876	H=1.E0/((double)LTF);
1877	naf_zardyd(ZTF,LTF,H,&ZA);
1878	*A = ZA.reel;
1879	*B = ZA.imag;
1880	*RM=module(ZA);
1881	naf_ztder(g_NAFVariable.KTABS,NVECT,ZTF,g_NAFVariable.ZTABS,TWIN,ZINC,ZEX,g_NAFVariable.T0,g_NAFVariable.XH); /CALL ZTDER(g_NAFVariable.KTABS+1,NVECT,ZTF,g_NAFVariable.ZTABS,TWIN,ZINC,ZEX,g_NAFVariable.T0,g_NAFVariable.XH)/
1882	naf_zardyd(ZTF,LTF,H,&ZB);
1883	DER=(mulcomplexe(conjcomplexe(ZA),ZB)).imag2.E0;
1884	#else /remplacee par: /
1885	/*!
1886	!------------------ CONVERSION DE FS EN RD/AN*/
1887	OM=FS/g_NAFVariable.UNIANG;
1888	i_compl_cmplx(&ZI,0.E0,1.E0);
1889	/*!------------------ LONGUEUR DU TABLEAU A INTEGRER
1890	!------------------ DANS LE CAS REEL MEME CHOSE */
1891	LTF=g_NAFVariable.KTABS;
1892	ANG0=OM*g_NAFVariable.T0;
1893	ANGI=OM*g_NAFVariable.XH;
1894	ZAC = i_compl_exp (i_compl_muldoubl(-ANG0,ZI));
1895	ZINC= i_compl_exp (i_compl_muldoubl(-ANGI,ZI));
1896	ZEX = i_compl_div(ZAC,ZINC); /ZEX = ZAC/ZINC/
1897	naf_ztpow2(g_NAFVariable.KTABS,NVECT,ZTF,g_NAFVariable.ZTABS,TWIN,ZINC,ZEX); /CALL ZTPOW2(g_NAFVariable.KTABS+1,NVECT,ZTF,g_NAFVariable.ZTABS,TWIN,ZINC,ZEX)/
1898	/!------------------ TAILLE DU PAS/
1899	H=1.E0/((double)LTF);
1900	naf_zardyd(ZTF,LTF,H,&ZA);
1901	*A = ZA.reel;
1902	*B = ZA.imag;
1903	*RM=i_compl_module(ZA);
1904	naf_ztder(g_NAFVariable.KTABS,NVECT,ZTF,g_NAFVariable.ZTABS,TWIN,ZINC,ZEX,g_NAFVariable.T0,g_NAFVariable.XH); /CALL ZTDER(g_NAFVariable.KTABS+1,NVECT,ZTF,g_NAFVariable.ZTABS,TWIN,ZINC,ZEX,g_NAFVariable.T0,g_NAFVariable.XH)/
1905	naf_zardyd(ZTF,LTF,H,&ZB);
1906	DER=(i_compl_mul(i_compl_conj(&ZA),ZB)).imag2.E0;
1907	#endif /NAF_USE_OPTIMIZE==0/
1908	/* v0.96 M. GASTINEAU 01/12/98 : fin optimisation */
1909	SYSFREE(ZTF);
1910	#undef NVECT
1911	}/* END SUBROUTINE PRODER*/
1912
1913	/* SUBROUTINE ZTDER (N,N1,ZTF,ZTA,TW,ZA,ZAST,T0,XH)*/
1914	/v0.96 M. GASTINEAU 09/09/98 : modification dans les boucles de I en I-1 /
1915	void naf_ztder(int N, int N1, t_complexe ZTF, t_complexe ZTA, double *TW, t_complexe ZA, t_complexe ZAST, double T0, double XH)
1916	{
1917	/*!-----------------------------------------------------------------------
1918	!ZTPOW CALCULE ZTF(I) = ZTA(I)* TW(I)ZASTZA*I (T0+(I-1)*XH EN VECTORIEL
1919	!
1920	!-----------------------------------------------------------------------
1921	IMPLICIT NONE
1922	! (N,N1,ZT,ZTF,ZTA,TW,ZA,ZAST,T0,XH)
1923	integer :: N,N1
1924	complex (8) :: ZTF(0:N),ZTA(0:N),ZA,ZAST
1925	REAL (8) :: TW(0:N),T0,XH
1926	!
1927	integer :: I,INC,NT,IT,NX
1928	complex (8) :: ZT1,ZINC
1929	complex (8),dimension (:), allocatable :: ZT
1930	! */
1931	int I,INC,NT,IT,NX;
1932	t_complexe ZT1, ZINC;
1933	t_complexe *ZT=NULL;
1934
1935	SYSCHECKMALLOCSIZE(ZT,t_complexe,N1); /allocate(ZT(0:N1-1),stat = NERROR)/
1936	if (N<N1-1)
1937	{
1938	fprintf(stdout,"DANS ZTDER, N = %d\n", N);
1939	return;
1940	}
1941	/!----------! /
1942	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
1943	#if NAF_USE_OPTIMIZE==0
1944	ZT[0] = mulcomplexe(ZAST,ZA);
1945	for( I = 1; I<N1; I++)
1946	{
1947	ZT[I] = mulcomplexe(ZT[I-1],ZA);
1948	}
1949	for(I = 0; I<N1; I++)
1950	{
1951	/v0.96M. GASTINEAU 09/09/98 : modification /
1952	/ZTF[I] = muldoublcomplexe((T0+(I-1)XH)TW[I], mulcomplexe(ZTA[I],ZT[I]));/ /ZTF(I) = ZTA(I)TW(I)ZT(I)(T0+(I-1)XH)/
1953	ZTF[I] = muldoublcomplexe((T0+(I)XH)TW[I], mulcomplexe(ZTA[I],ZT[I])); /ZTF(I) = ZTA(I)TW(I)ZT(I)(T0+(I-1)XH)/
1954	}
1955	ZT1 = divcomplexe(ZT[N1-1],ZAST);
1956	ZINC= cmplx(1,0);
1957	INC =0;
1958	NT = (N+1)/N1;
1959	for( IT = 2; IT<= NT; IT++)
1960	{
1961	ZINC = mulcomplexe(ZINC,ZT1);
1962	INC += N1; /INC = INC + N1/
1963	for(I = 0; I<N1; I++)
1964	{
1965	/v0.96M. GASTINEAU 09/09/98 : modification /
1966	/ZTF[INC+I]=muldoublcomplexe((T0+(INC+I-1)XH)TW[INC+I],mulcomplexe(ZTA[INC+I],mulcomplexe(ZT[I],ZINC)));/
1967	ZTF[INC+I]=muldoublcomplexe((T0+(INC+I)XH)TW[INC+I],mulcomplexe(ZTA[INC+I],mulcomplexe(ZT[I],ZINC)));
1968	/ZTF(INC+I)=ZTA(INC+I)TW(INC+I)ZT(I)ZINC(T0+(INC+I-1)XH)*/
1969	}
1970	}
1971	ZINC = mulcomplexe(ZINC,ZT1);
1972	INC += N1; /INC = INC + N1/
1973	NX = (N+1)-NT*N1;
1974	for (I = 0; I<NX; I++)
1975	{
1976	/v0.96M. GASTINEAU 09/09/98 : modification /
1977	/ZTF[INC+I]=muldoublcomplexe(TW[INC+I](T0+(INC+I-1)XH), mulcomplexe(ZTA[INC+I], mulcomplexe(ZT[I],ZINC)));/
1978	ZTF[INC+I]=muldoublcomplexe(TW[INC+I](T0+(INC+I)XH), mulcomplexe(ZTA[INC+I], mulcomplexe(ZT[I],ZINC)));
1979	/ZTF(INC+I)=ZTA(INC+I)TW(INC+I)ZT(I)ZINC(T0+(INC+I-1)XH)*/
1980	}
1981	#else /remplacee par:/
1982	ZT[0] = i_compl_mul(ZAST,ZA);
1983	for( I = 1; I<N1; I++)
1984	{
1985	ZT[I] = i_compl_mul(ZT[I-1],ZA);
1986	}
1987	for(I = 0; I<N1; I++)
1988	{
1989	ZTF[I] = i_compl_muldoubl((T0+(I)XH)TW[I], i_compl_mul(ZTA[I],ZT[I]));
1990	/ZTF(I) = ZTA(I)TW(I)ZT(I)(T0+(I-1)XH)/
1991	}
1992	ZT1 = i_compl_div(ZT[N1-1],ZAST);
1993	i_compl_cmplx(&ZINC,1,0);
1994	INC =0;
1995	NT = (N+1)/N1;
1996	for( IT = 2; IT<= NT; IT++)
1997	{
1998	i_compl_pmul(&ZINC,&ZT1);
1999	INC += N1; /INC = INC + N1/
2000	for(I = 0; I<N1; I++)
2001	{
2002	ZTF[INC+I]=i_compl_muldoubl((T0+(INC+I)XH)TW[INC+I],i_compl_mul(ZTA[INC+I],i_compl_mul(ZT[I],ZINC)));
2003	/ZTF(INC+I)=ZTA(INC+I)TW(INC+I)ZT(I)ZINC(T0+(INC+I-1)XH)*/
2004	}
2005	}
2006	i_compl_pmul(&ZINC,&ZT1);
2007	INC += N1; /INC = INC + N1/
2008	NX = (N+1)-NT*N1;
2009	for (I = 0; I<NX; I++)
2010	{
2011	ZTF[INC+I]=i_compl_muldoubl(TW[INC+I](T0+(INC+I)XH), i_compl_mul(ZTA[INC+I], i_compl_mul(ZT[I],ZINC)));
2012	/ZTF(INC+I)=ZTA(INC+I)TW(INC+I)ZT(I)ZINC(T0+(INC+I-1)XH)*/
2013	}
2014	#endif /NAF_USE_OPTIMIZE==0/
2015	/* v0.96 M. GASTINEAU 01/12/98 : fin optimisation */
2016	SYSFREE(ZT);
2017	}/* END SUBROUTINE ZTDER*/
2018
2019	/* SUBROUTINE SECANTES(X,PASS,EPS,XM,IPRT,NFPRT)*/
2020	void naf_secantes(double X, double PASS, double EPS, double XM, int IPRT, FILE NFPRT)
2021	{
2022	/*! SUBROUTINE SECANTES(X,PASS,EPS,XM,FONC,IPRT,NFPRT)
2023	!***********************************************************
2024	! CALCUL PAR LA METHODE DES SECANTES D'UN ZERO
2025	! D'UNE FONCTION FUNC FOURNIE PAR L'UTILISATEUR
2026	!
2027	! X, PASS : LE ZERO EST SUPPOSE ETRE DANS X-PASS,X+PASS
2028	! EPS : PRECISION AVEC LAQUELLE ON VA CALCULER
2029	! LE MAXIMUM
2030	! XM : ABSCISSE DU ZERO
2031	! g_NAFVariable.IPRT : -1 PAS D'IMPRESSION
2032	! 0: IMPRESSION D'UNE EVENTUELLE ERREUR FATALE
2033	! 1: IMPRESSION SUR LE FICHIER g_NAFVariable.NFPRT
2034	! 2: IMPRESSION DES RESULTATS INTERMEDIAIRES
2035	!
2036	! FONC(X) :FONCTION DONNEE PAR L'UTILISATEUR
2037	! PARAMETRES
2038	! NMAX : NOMBRE D'ESSAIS
2039	! IENC : 1 : test d'encadrement 0 sinon
2040	! VARIABLE NERROR
2041	! 0: A PRIORI TOUT S'EST BIEN PASSE
2042	! 1: LE ZERO N'EST PAS TROUVE A EPSI PRES
2043	! --> SEULEMENT SI IENC=1
2044	! 2: PENTRE TROP FAIBLE (DIVISION PAR PRESQUE ZERO)
2045	! 3: ECHEC TOTAL
2046	! F. Joutel 26/5/98
2047	! Modif 26/8/98 Enleve l'appel de FONC car dans un module
2048	! il ne semble pas etre possible de passer une fonction
2049	! interne en parametre, remplace FONC par FUNCP
2050	!***********************************************************
2051
2052	IMPLICIT NONE
2053	! (X,PASS,EPS,XM,FUNCP,IPRT,NFPRT)
2054	real (8) :: X, PASS,EPS,XM
2055	integer :: IPRT,NFPRT
2056	! interface
2057	! function fonc(x)
2058	! real (8) :: x,fonc
2059	! end function fonc
2060	! end interface
2061	!
2062	integer, parameter :: NMAX=30, IENC=1
2063	integer :: I
2064	REAL (8) :: EPSI,A,B,FA,FB,DELTA,COR
2065	!*/
2066	#define SECANTES_NMAX (30)
2067	#define SECANTES_IENC (1)
2068	int I;
2069	double EPSI,A,B,FA,FB,DELTA,COR;
2070
2071	g_NAFVariable.NERROR=0;
2072	EPSI=MAX(g_NAFVariable.EPSM,EPS);
2073	if (fabs(PASS)>EPSI)
2074	{
2075	if (IPRT>=1)
2076	{
2077	fprintf(NFPRT," AMELIORATION PAR LES SECANTES\n");
2078	}
2079	}
2080	I=0;
2081	A=X-PASS;
2082	B=X;
2083	FB=naf_funcp(A);
2084	/! FB=FONC(A)/
2085	/10 CONTINUE/
2086	SECANTES_10 :
2087	if (fabs(B-A)>EPSI)
2088	{
2089	FA=FB;
2090	FB=naf_funcp(B);
2091	/! FB=FONC(B)/
2092	DELTA= FB-FA;
2093	if (IPRT>=2)
2094	{
2095	fprintf(NFPRT,"SEC: A=%g, B=%g, abs(B-A)=%g \n", A, B,fabs(B-A));
2096	fprintf(NFPRT,"SEC: F(A)=%g, F(B)=%g\n", FA, FB);
2097	}
2098	if (fabs(DELTA)<=g_NAFVariable.EPSM)
2099	{
2100	g_NAFVariable.NERROR=2;
2101	if (IPRT>=1)
2102	{
2103	fprintf(NFPRT,"ECHEC DE LA METHODE DES SECANTES\n");
2104	fprintf(NFPRT," DIVISION PAR PRESQUEZERO\n");
2105	fprintf(NFPRT," ON CONTINUE AVEC LA VALEUR TROUVEE\n");
2106	}
2107	*XM=B;
2108	return;
2109	}
2110	COR = FB*(A-B)/DELTA;
2111	A = B;
2112	B = B+COR;
2113	I++; /I=I+1 /
2114	if (I>SECANTES_NMAX)
2115	{
2116	g_NAFVariable.NERROR=3;
2117	if (IPRT>=0)
2118	{
2119	fprintf(NFPRT,"ECHEC DE LA METHODE DES SECANTES\n");
2120	fprintf(NFPRT," BEAUCOUP TROP D ITERATIONS\n");
2121	fprintf(NFPRT," ON CONTINUE AVEC LA VALEUR INITIALE\n");
2122	}
2123	*XM=X;
2124	return;
2125	}
2126	goto SECANTES_10; /GOTO 10/
2127	}
2128	/!----- !fin de la boucle/
2129	*XM=B;
2130	if (SECANTES_IENC==1)
2131	{
2132	if (naf_funcp(B-EPSI)*naf_funcp(B+EPSI)>0.E0)
2133	{
2134	/! if (FONC(B-EPSI)FONC(B+EPSI)>0.D0) {*/
2135	g_NAFVariable.NERROR=1;
2136	}
2137	}
2138	if (IPRT==1)
2139	{
2140	/v0.96 M. GASTINEAU 19/11/98 : modification/
2141	/fprintf(NFPRT,"POSITION SECANTES %g TROUVEE A %g",XM, EPSI);/
2142	/remplacee par:/
2143	fprintf(NFPRT,"POSITION SECANTES %g TROUVEE A %g",*XM, EPSI);
2144	/v0.96 M. GASTINEAU 19/11/98 : fin modification/
2145	if (SECANTES_IENC==1)
2146	{
2147	if (g_NAFVariable.NERROR==1)
2148	{
2149	fprintf(NFPRT," SANS GARANTIE\n");
2150	}
2151	if (g_NAFVariable.NERROR==0)
2152	{
2153	fprintf(NFPRT," AVEC GARANTIE\n");
2154	}
2155	}
2156	/v0.96 M. GASTINEAU 19/11/98 : modification/
2157	/* fprintf(NFPRT,"\n %19.9E %12.8E\n", PASS,XM);*/
2158	/remplacee par:/
2159	fprintf(NFPRT,"\n %19.9E %12.8E\n", PASS,*XM);
2160	/v0.96 M. GASTINEAU 19/11/98 : fin modification/
2161	}
2162	#undef SECANTES_NMAX
2163	#undef SECANTES_IENC
2164	}/* END SUBROUTINE SECANTES*/
2165
2166	/* SUBROUTINE MAXIQUA(X,PASS,EPS,XM,YM,IPRT,NFPRT)*/
2167	void naf_maxiqua(double X, double PASS, double EPS, double XM, double YM, int IPRT, FILE *NFPRT)
2168	{
2169	/*
2170	! SUBROUTINE MAXIQUA(X,PASS,EPS,XM,YM,FONC,IPRT,NFPRT)
2171	!***********************************************************
2172	! CALCUL PAR INTERPOLATION QUADRATIQUE DU MAX
2173	! D'UNE FONCTION FONC FOURNIE PAR L'UTILISATEUR
2174	!
2175	! X, PASS : LE MAX EST SUPPOSE ETRE SI POSSIBLE
2176	! DANS X-PASS,X+PASS
2177	! EPS : PRECISION AVEC LAQUELLE ON VA CALCULER
2178	! LE MAXIMUM
2179	! XM : ABSCISSE DU MAX
2180	! YM : YM=FUNC(XM)
2181	! IPRT : -1 PAS D'IMPRESSION
2182	! 0: IMPRESSION D'UNE EVENTUELLE ERREUR FATALE
2183	! 1: IMPRESSION SUR LE FICHIER NFPRT
2184	! 2: IMPRESSION DES RESULTATS INTERMEDIAIRES
2185	!
2186	! FONC(X) :FONCTION DONNEE PAR L'UTILISATEUR
2187	!
2188	!***********************************************************
2189	! PARAMETRES
2190	! NMAX : NOMBRE MAXIMAL D'ITERATIONS
2191	! VARIABLE (COMMON ASD)
2192	! g_NAFVariable.NERROR:
2193	! 0 : OK
2194	! 1 : COURBURE BIEN FAIBLE !
2195	! 2 : LA FONCTION EST BIEN PLATE !
2196	! 3 : OSCILLATION DE L'ENCADREMENT
2197	! 4 : PAS D'ENCADREMENT TROUVE
2198	! 5 : ERREUR FATALE
2199	!***********************************************************
2200	! ASD VERSION 0.9 J.LASKAR & F.JOUTEL 16/5/97
2201	! CHANGEMENT POUR UN EPSILON MACHINE CALCULE PAR CALCEPSM
2202	! JXL 22/4/98
2203	! Modif 26/8/98 Enleve l'appel de FONC car dans un module
2204	! il ne semble pas etre possible de passer une fonction
2205	! interne en parametre, remplace FONC par FUNCP
2206	!***********************************************************
2207
2208	IMPLICIT NONE
2209	! (X,PASS,EPS,XM,YM,FONC,IPRT,NFPRT)
2210	integer :: IPRT,NFPRT
2211	REAL (8) :: X,PASS,EPS,XM,YM
2212
2213	! interface
2214	! function fonc (x)
2215	! real (8) :: x,fonc
2216	! end function fonc
2217	! end interface
2218	!
2219	integer, parameter :: NMAX=200, NFATAL=100
2220	REAL (8), parameter :: RABS=1.D0
2221	integer :: NITER,NCALCUL
2222	REAL (8) :: PAS,EPSLOC,X1,X2,X3,Y1,Y2,Y3,A,ERR,DX,TEMP
2223	REAL (8) :: D1,D2
2224	*/
2225	#define MAXIQUA_NMAX (200)
2226	#define MAXIQUA_NFATAL (100)
2227	#define MAXIQUA_RABS (1.E0)
2228	int NITER, NCALCUL;
2229	double PAS,EPSLOC,X1,X2,X3,Y1,Y2,Y3,A,ERR,DX,TEMP;
2230	double D1=0E0,D2=0E0;
2231
2232	/*!
2233	!* ENTRE 0 ET RABS LES ERREURS SONT ENVISAGEES DE FACON ABSOLUE
2234	!* APRES RABS, LES ERREURS SONT ENVISAGEES DE FACON RELATIVE.
2235	!*/
2236	EPSLOC=MAX(EPS,sqrt(g_NAFVariable.EPSM));
2237	/*!
2238	!* AU SOMMET D'UNE PARABOLE Y=Y0-A(X-X0)*2, UN ECART DE X AUTOUR DE X0
2239	!* INFERIEUR A EPSLOC DONNE UN ECART DE Y INFERIEUR A AEPSLOC2/
2240	if (IPRT>=1)
2241	{
2242	fprintf(NFPRT," ROUTINE DE RECHERCHE DU MAXIMUM :\n");
2243	}
2244	PAS = PASS;
2245	NITER =0;
2246	NCALCUL=0;
2247	g_NAFVariable.NERROR=0;
2248	X2 = X;
2249	Y2 = naf_func (X2);
2250	X1 = X2 - PAS;
2251	X3 = X2 + PAS;
2252	Y1 = naf_func (X1);
2253	Y3 = naf_func (X3);
2254	/! DECALAGE POUR ENCADRER LE MAXIMUM*/
2255	MAXIQUA_10:
2256	if ((NITER>MAXIQUA_NMAX)\|\|(NCALCUL>MAXIQUA_NFATAL))
2257	{
2258	if (NCALCUL>MAXIQUA_NFATAL)
2259	{
2260	g_NAFVariable.NERROR =5;
2261	if (IPRT>=0)
2262	{
2263	fprintf(NFPRT," ERREUR FATALE\n");
2264	}
2265	}
2266	else
2267	{
2268	if (PAS>=PASS)
2269	{
2270	g_NAFVariable.NERROR = 4;
2271	if (IPRT>=0)
2272	{
2273	fprintf(NFPRT," PAS D''ENCADREMENT TROUVE\n");
2274	}
2275	}
2276	else
2277	{
2278	g_NAFVariable.NERROR =3;
2279	if (IPRT>=0)
2280	{
2281	fprintf(NFPRT," OSCILLATION DE L''ENCADREMENT ?\n");
2282	}
2283	}
2284	}
2285	ERR=PAS;
2286	goto MAXIQUA_EXIT;
2287	}
2288	if ((Y1>Y2)\|\|(Y3>Y2))
2289	{
2290	NITER ++; /NITER = NITER +1/
2291	if (Y1>Y3)
2292	{
2293	X2 = X1;
2294	Y2 = Y1;
2295	}
2296	else
2297	{
2298	X2 = X3;
2299	Y2 = Y3;
2300	}
2301	X1 = X2 - PAS;
2302	X3 = X2 + PAS;
2303	Y1 = naf_func (X1);
2304	Y3 = naf_func (X3);
2305	goto MAXIQUA_10;
2306	}
2307	D1 = Y2-Y1;
2308	D2 = Y3-Y2;
2309	A = fabs(0.5E0(D2-D1)/(PASPAS));
2310	if (IPRT>=2)
2311	{
2312	if (sqrt(A)!=0)
2313	{
2314	fprintf(NFPRT,"%12.8E %12.8E %25.16E %25.16E %25.16E %25.16E\n", PAS,X2,Y1,Y2,Y3,EPSLOC/sqrt(A));
2315	}
2316	else
2317	{
2318	fprintf(NFPRT,"%12.8E %12.8E %25.16E %25.16E %25.16E INF\n", PAS,X2,Y1,Y2,Y3);
2319	}
2320	}
2321	/! TEST POUR UN CALCUL SIGNifICATif DES VALEURS DE LA FONCTION*/
2322	if ((fabs (D1-D2))<2.E0g_NAFVariable.EPSMMAX(MAXIQUA_RABS,fabs(Y2)))
2323	{
2324	g_NAFVariable.NERROR= 2;
2325	if (IPRT>=1)
2326	{
2327	fprintf(NFPRT," PLATEAU DE LA FONCTION\n");
2328	}
2329	ERR=PAS;
2330	goto MAXIQUA_EXIT;
2331	}
2332
2333	/! TEST SUR LA FORME DE LA COURBE */
2334	if (A<g_NAFVariable.EPSM*MAX(MAXIQUA_RABS,fabs(Y2)))
2335	{
2336	g_NAFVariable.NERROR= 1;
2337	if (IPRT>=1)
2338	{
2339	fprintf(NFPRT,"COURBE TROP PLATE\n");
2340	}
2341	ERR=PAS;
2342	goto MAXIQUA_EXIT;
2343	}
2344	/! TEST POUR UN CALCUL SIGNifICATif AU SOMMET D'UNE PARABOLE
2345	!* ON PEUT IGNORER LE TEST MAIS ON N'A PLUS DE GARANTIE SUR
2346	!* LES VALEURS DE LA FONCTION PLUS TARD */
2347	if (PAS<=EPSLOC*sqrt(MAX(MAXIQUA_RABS,fabs(Y2))/A))
2348	{
2349	if (IPRT>=1)
2350	{
2351	fprintf(NFPRT,"SORTIE AU SOMMET DE LA PARABOLE\n");
2352	}
2353	ERR=EPSLOC*sqrt(MAX(MAXIQUA_RABS,Y2)/A);
2354	goto MAXIQUA_EXIT;
2355	}
2356	DX= 0.5E0PAS(Y1-Y3)/(D2-D1);
2357	/! TEST POUR UN NOUVEAU PAS SIGNifICATif*/
2358	if (fabs(DX)<g_NAFVariable.EPSM*MAX(MAXIQUA_RABS,fabs(X2)))
2359	{
2360	if (IPRT>=1)
2361	{
2362	fprintf(NFPRT,"CORRECTION MINUSCULE\n");
2363	}
2364	ERR=PAS;
2365	goto MAXIQUA_EXIT;
2366	}
2367	PAS=fabs(DX);
2368	if (DX>0.E0)
2369	{
2370	X1 = X2;
2371	Y1 = Y2;
2372	X2 +=PAS; /X2 = X2+PAS/
2373	Y2 = naf_func(X2);
2374	X3 = X2 + PAS;
2375	Y3 = naf_func(X3);
2376	}
2377	else
2378	{
2379	X3 = X2;
2380	Y3 = Y2;
2381	X2 -=PAS; /X2 = X2-PAS/
2382	Y2 = naf_func(X2);
2383	X1 = X2 - PAS;
2384	Y1 = naf_func(X1);
2385	}
2386	/! A LA SORTIE DE CE CALCUL LE PAS DOIT ETRE DIVISE AU MOINS PAR 2
2387	!! MAIS ON N'EST PAS ASSURE QUE LES VALEURS FINALES Y1,Y2,Y3
2388	!* VERifIENT Y1<=Y2>=Y3, AUSSI, ON REPART AU DEBUT*/
2389	NCALCUL++; /NCALCUL=NCALCUL+1/
2390	goto MAXIQUA_10;
2391	MAXIQUA_EXIT :
2392	*XM = X2;
2393	*YM = Y2;
2394	if (IPRT>=1)
2395	{
2396	fprintf(NFPRT,"%12.8E %12.8E %12.8E\n", PAS,XM,YM);
2397	fprintf(NFPRT," POSITION TROUVEE A %g PRES\n",ERR);
2398	TEMP = fabs (Y3-Y2)+ fabs(Y2-Y1);
2399	if (TEMP<2.E0g_NAFVariable.EPSMMAX(MAXIQUA_RABS,fabs(Y2)))
2400	{
2401	fprintf(NFPRT," TROUVEE SOUS UN PLATEAU DE LA FONCTION\n");
2402	}
2403	}
2404	#undef MAXIQUA_NMAX
2405	#undef MAXIQUA_NFATAL
2406	#undef MAXIQUA_RABS
2407	}/* END SUBROUTINE MAXIQUA*/
2408
2409	/* SUBROUTINE FREFIN (FR,A,B,RM, RPAS0,RPREC)*/
2410	/v0.96 M. GASTINEAU 07/12/98 : modification car bug lors de l'optimisation sur Mac /
2411	void naf_frefin(double FR, double A, double B, double RM, const double RPAS0, const double RPREC)
2412	{
2413	/*!-----------------------------------------------------------------------
2414	! FREFIN RECHERCHE FINE DE LA FREQUENCE
2415	! FR FREQUENCE DE BASE EN "/AN ENTREE ET SORT
2416	! A PARTIE REELLE DE L'AMPLITUDE (SORTIE)
2417	! B PARTIE IMAGINAIRE DE L'AMPLITUDE (SORTIE)
2418	! RM MODULE DE L'AMPLITUDE (SORTIE)
2419	! RPAS0 PAS INITIAL EN "/AN
2420	! RPREC PRECISION DEMANDEE EN "/AN
2421	! REVU LE 26/9/87 J. LASKAR
2422	! MODif LE 15/11/90 MAX PAR INTERP. QUAD.
2423	! Modif 26/8/98 Enleve les appels aux fonctions func et funcp internes
2424	! au module
2425	!-----------------------------------------------------------------------
2426
2427	IMPLICIT NONE
2428	! (FR,A,B,RM, RPAS0,RPREC)
2429	REAL (8) FR,A,B,RM,RPAS0,RPREC
2430	!
2431	REAL (8) X,PASS,EPS,XM,YM
2432	! */
2433	double X,PASS,EPS,XM,YM;
2434	/v0.96 M. GASTINEAU 07/12/98 : modification car bug quand optimisation sur Mac /
2435	/* X = *FR;
2436	PASS = RPAS0;
2437	EPS = RPREC;/ /remplacee par:*/
2438	PASS = RPAS0;
2439	EPS = RPREC;
2440	X = *FR;
2441	/v0.96 M. GASTINEAU 07/12/98 : fin modification /
2442	naf_maxiqua(X,PASS,EPS,&XM,&YM,g_NAFVariable.IPRT,g_NAFVariable.NFPRT);
2443	/*! CALL MAXIQUA(X,PASS,EPS,XM,YM,FUNC,g_NAFVariable.IPRT,g_NAFVariable.NFPRT)
2444	!************! AMELIORATION DE LA RECHERCHE PAR LE METHODE DES SECANTES
2445	!***********! APPLIQUEE A LA DERIVEE DU MODULE./
2446	if (g_NAFVariable.ISEC==1)
2447	{
2448	X=XM;
2449	naf_secantes(X,PASS,EPS,&XM,g_NAFVariable.IPRT,g_NAFVariable.NFPRT);
2450	/! CALL SECANTES(X,PASS,EPS,XM,FUNCP,g_NAFVariable.IPRT,g_NAFVariable.NFPRT)/
2451	YM=naf_func(XM);
2452	}
2453	*FR=XM;
2454	*RM=YM;
2455	*A=AF;
2456	*B=BF;
2457	if (g_NAFVariable.IPRT==1)
2458	{
2459	/v0.96 M. GASTINEAU 19/11/98 : modification/
2460	/fprintf(g_NAFVariable.NFPRT,"\n%10.6E %19.9E %19.9E %19.9E\n", FR,RM,AF,BF);/
2461	/remplacee par:/
2462	fprintf(g_NAFVariable.NFPRT,"\n%10.6E %19.9E %19.9E %19.9E\n", FR,RM,AF,BF);
2463	/v0.96 M. GASTINEAU 19/11/98 : fin modification/
2464	}
2465	/1000 FORMAT (1X,F10.6,3D19.9)/
2466	}/* END SUBROUTINE FREFIN*/
2467
2468
2469
2470	/* SUBROUTINE PROFRE(FS,A,B,RMD)*/
2471	BOOL naf_profre(double FS, double A, double B, double *RMD)
2472	{
2473	/!**********************************************************************
2474	! CE SOUS PROG. CALCULE LE PRODUIT SCALAIRE DE EXP(-IFS)T PAR TAB(0:7
2475	! UTILISE LA METHODE D'INTEGRATION DE HARDY
2476	! LES RESULTATS SONT DANS A,B,RMD PARTIE REELLE, IMAGINAIRE, MODULE
2477	! FS EST DONNEE EN SECONDES PAR AN
2478	! REVU LE 26/9/87 J. LASKAR
2479	!***********************************************************************
2480
2481	IMPLICIT NONE
2482	! (g_NAFVariable.KTABS,FS,A,B,RMD)
2483	REAL (8) :: FS,A,B,RMD
2484	!
2485	integer :: LTF
2486	REAL (8) :: OM,ANG0,ANGI,H
2487	complex (8) ZI,ZAC,ZINC,ZEX,ZA
2488	complex (8), dimension (:), allocatable :: ZTF
2489	! */
2490	int LTF;
2491	double OM,ANG0,ANGI,H;
2492	t_complexe ZI,ZAC,ZINC,ZEX,ZA;
2493	t_complexe *ZTF=NULL;
2494
2495	SYSCHECKMALLOCSIZE(ZTF, t_complexe, g_NAFVariable.KTABS+1); /allocate(ZTF(0:g_NAFVariable.KTABS),stat = g_NAFVariable.NERROR)/
2496	/*!
2497	!------------------ CONVERSION DE FS EN RD/AN*/
2498	OM=FS/g_NAFVariable.UNIANG ;
2499	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
2500	#if NAF_USE_OPTIMIZE==0
2501	ZI=cmplx(0.E0,1.E0);
2502	#else /remplacee par:/
2503	i_compl_cmplx(&ZI,0.E0,1.E0);
2504	#endif /NAF_USE_OPTIMIZE==0/
2505	/* v0.96 M. GASTINEAU 01/12/98 : fin optimisation */
2506	/*!------------------ LONGUEUR DU TABLEAU A INTEGRER
2507	!------------------ DANS LE CAS REEL MEME CHOSE */
2508	LTF=g_NAFVariable.KTABS;
2509	ANG0=OM*g_NAFVariable.T0;
2510	ANGI=OM*g_NAFVariable.XH;
2511	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
2512	#if NAF_USE_OPTIMIZE==0
2513	ZAC = expcomplexe (muldoublcomplexe(-ANG0,ZI));
2514	ZINC= expcomplexe (muldoublcomplexe(-ANGI,ZI));
2515	ZEX = divcomplexe(ZAC,ZINC);
2516	#else /remplacee par:/
2517	ZAC = i_compl_exp(i_compl_muldoubl(-ANG0,ZI));
2518	ZINC= i_compl_exp(i_compl_muldoubl(-ANGI,ZI));
2519	ZEX = i_compl_div(ZAC,ZINC);
2520	#endif /NAF_USE_OPTIMIZE==0/
2521
2522	/* v0.96 M. GASTINEAU 01/12/98 : fin optimisation */
2523	naf_ztpow2(g_NAFVariable.KTABS,64,ZTF,g_NAFVariable.ZTABS,TWIN,ZINC,ZEX); /CALL ZTPOW2(g_NAFVariable.KTABS+1,64,ZTF,g_NAFVariable.ZTABS,TWIN,ZINC,ZEX)/
2524	/!------------------ TAILLE DU PAS/
2525	H=1.E0/((double)LTF);
2526	if (naf_zardyd(ZTF,LTF,H,&ZA)==FALSE)
2527	{
2528	SYSFREE(ZTF);
2529	return FALSE;
2530	}
2531	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
2532	#if NAF_USE_OPTIMIZE==0
2533	*RMD=module(ZA);
2534	#else /remplacee par:/
2535	*RMD=i_compl_module(ZA);
2536	#endif /NAF_USE_OPTIMIZE==0/
2537	/* v0.96 M. GASTINEAU 01/12/98 : fin optimisation */
2538	*A = ZA.reel;
2539	*B = ZA.imag;
2540	SYSFREE(ZTF);
2541	return TRUE;
2542	}/* END SUBROUTINE PROFRE*/
2543
2544	/* SUBROUTINE ZTPOW2 (N,N1,ZTF,ZTA,TW,ZA,ZAST)*/
2545	void naf_ztpow2(int N, int N1, t_complexe ZTF, t_complexe ZTA, double *TW, t_complexe ZA, t_complexe ZAST)
2546	{
2547	/*!-----------------------------------------------------------------------
2548	! ZTPOW CALCULE ZTF(I) = ZTA(I)* TW(I)ZASTZA**I EN VECTORIEL
2549	!
2550	!-----------------------------------------------------------------------
2551	IMPLICIT NONE
2552	! (N,N1,ZTF,ZTA,TW,ZA,ZAST)
2553	integer :: N,N1
2554	complex (8) :: ZTF(0:N),ZTA(0:N),ZA,ZAST
2555	REAL (8) :: TW(0:N)
2556	!
2557	integer :: INC,NT,IT,I,NX
2558	complex (8) :: ZT1,ZINC
2559	complex (8), dimension(:),allocatable :: ZT
2560	! */
2561	int INC,NT,IT,I,NX;
2562	t_complexe ZT1,ZINC;
2563	t_complexe *ZT=NULL;
2564	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
2565	#if NAF_USE_OPTIMIZE>0
2566	t_complexe *pzarZT;
2567	#endif /**/
2568	SYSCHECKMALLOCSIZE(ZT,t_complexe, N1); /allocate(ZT(0:N1-1),stat = nerror)/
2569	if (N<N1-1)
2570	{
2571	fprintf(stdout,"DANS ZTPOW, N = %d\n", N);
2572	return;
2573	}
2574	/!----------! /
2575	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
2576	#if NAF_USE_OPTIMIZE==0
2577	ZT[0] = mulcomplexe(ZAST,ZA);
2578	for(I = 1; I<N1; I++)
2579	{
2580	ZT[I]=mulcomplexe(ZT[I-1], ZA); /ZT(I) = ZT(I-1)ZA*/
2581	}
2582	for(I = 0; I<N1; I++)
2583	{
2584	ZTF[I]=muldoublcomplexe(TW[I], mulcomplexe(ZTA[I], ZT[I])); /ZTF(I) = ZTA(I)TW(I)ZT(I)/
2585	}
2586	ZT1= divcomplexe(ZT[N1-1], ZAST); /ZT1 = ZT(N1-1)/ZAST/
2587	ZINC= cmplx(1E0,0E0);
2588	INC =0;
2589	NT = (N+1)/N1;
2590	for(IT = 2;IT<= NT; IT++)
2591	{
2592	ZINC=mulcomplexe(ZINC,ZT1); /ZINC = ZINCZT1*/
2593	INC += N1; /INC = INC + N1/
2594	for(I = 0; I< N1;I++)
2595	{
2596	/ZTF(INC +I) = ZTA(INC+I)TW(INC+I)ZT(I)ZINC*/
2597	ZTF[INC +I] = muldoublcomplexe(TW[INC+I], mulcomplexe(mulcomplexe(ZTA[INC+I], ZT[I]),ZINC));
2598	}
2599	}
2600	ZINC = mulcomplexe(ZINC, ZT1); /ZINC = ZINCZT1*/
2601	INC += N1; /INC = INC + N1/
2602	NX = N+1-NT*N1;
2603	for(I = 0; I< NX; I++)
2604	{
2605	/* ZTF(INC +I) = ZTA(INC+I)TW(INC+I)ZT(I)ZINC/
2606	ZTF[INC +I] = muldoublcomplexe(TW[INC+I], mulcomplexe(mulcomplexe(ZTA[INC+I], ZT[I]),ZINC));
2607	}
2608	#else /remplacee par:/
2609	ZT[0] = i_compl_mul(ZAST,ZA);
2610	for(I = 1, pzarZT=ZT+1; I<N1; I++, pzarZT++)
2611	{
2612	pzarZT=i_compl_mul((pzarZT-1), ZA); /ZT(I) = ZT(I-1)ZA*/
2613	}
2614	for(I = 0; I<N1; I++)
2615	{
2616	ZTF[I]=i_compl_muldoubl(TW[I], i_compl_mul(ZTA[I], ZT[I])); /ZTF(I) = ZTA(I)TW(I)ZT(I)/
2617	}
2618	ZT1= i_compl_div(ZT[N1-1], ZAST); /ZT1 = ZT(N1-1)/ZAST/
2619	i_compl_cmplx(&ZINC,1E0,0E0);
2620	INC =0;
2621	NT = (N+1)/N1;
2622	for(IT = 2;IT<= NT; IT++)
2623	{
2624	i_compl_pmul(&ZINC,&ZT1); /ZINC = ZINCZT1*/
2625	INC += N1; /INC = INC + N1/
2626	for(I = 0; I< N1;I++)
2627	{
2628	/ZTF(INC +I) = ZTA(INC+I)TW(INC+I)ZT(I)ZINC*/
2629	ZTF[INC +I] = i_compl_muldoubl(TW[INC+I], i_compl_mul(i_compl_mul(ZTA[INC+I], ZT[I]),ZINC));
2630	}
2631	}
2632	i_compl_pmul(&ZINC,&ZT1); /ZINC = ZINCZT1*/
2633	INC += N1; /INC = INC + N1/
2634	NX = N+1-NT*N1;
2635	for(I = 0; I< NX; I++)
2636	{
2637	/* ZTF(INC +I) = ZTA(INC+I)TW(INC+I)ZT(I)ZINC/
2638	ZTF[INC +I] = i_compl_muldoubl(TW[INC+I], i_compl_mul(i_compl_mul(ZTA[INC+I], ZT[I]),ZINC));
2639	}
2640	#endif /NAF_USE_OPTIMIZE==0/
2641	/* v0.96 M. GASTINEAU 01/12/98 : fin optimisation */
2642	SYSFREE(ZT);
2643	}/* END SUBROUTINE ZTPOW2*/
2644
2645	/* SUBROUTINE INIWIN*/
2646	/v0.96 M. GASTINEAU 18/12/98 : modification du prototype /
2647	/void naf_iniwin()//remplacee par: /
2648	void naf_iniwin(double *p_pardTWIN)
2649	/v0.96 M. GASTINEAU 18/12/98 : fin modification /
2650	{
2651	/!*****************************************************************
2652	! INITIALISE LE TABLEAU TWIN POUR UTILISATION D'UNE FENETRE
2653	! DANS L'ANALYSE DE FOURIER.
2654	! FENETRE DE HANNING:
2655	! PHI(T) = (1+COS(PI*T))
2656	! MODif LE 13/12/90 (J. LASKAR)
2657	! MODif LE 27/1/96 FOR VARIOUS WINDOWS (J. LASKAR)
2658	!
2659	! WORKING AREA TWIN(*) SHOULD BE GIVEN IN
2660	!
2661	! IW IS THE WINDOW FLAG
2662	! IW = 0 : NO WINDOW
2663	! IW = N > 0 PHI(T) = CN(1+COS(PIT))**N
2664	! WITH CN = 2^N*(N!)^2/(2N)!
2665	! IW = -1 EXPONENTIAL WINDOW PHI(T) = 1/CE*EXP(-1/(1-T^2))
2666	!
2667	! MODif LE 22/4/98 POUR CALCUL EN PLUS DE L'EPSILON MACHINE
2668	! g_NAFVariable.EPSM
2669	! 26/5/98 CORRECTION DE L'ORIGINE DES TABLEAUX (*m/4)
2670	!******************************************************************
2671	!
2672
2673	IMPLICIT NONE
2674	!
2675	integer :: I,IT
2676	REAL (8) :: CE,T1,T2,TM,T,CN,PIST
2677	!
2678	!---------------------------------------------------------------- */
2679	int I, IT;
2680	double CE,T1,T2,TM,T,CN,PIST;
2681
2682	CE= 0.22199690808403971891E0;
2683	/*!
2684	! PI=ATAN2(1.D0,0.E0)2.E0/
2685	T1=g_NAFVariable.T0;
2686	T2=g_NAFVariable.T0+g_NAFVariable.KTABS*g_NAFVariable.XH;
2687	TM=(T2-T1)/2;
2688	PIST=pi/TM;
2689	if (g_NAFVariable.IW==0)
2690	{
2691	for(IT=0;IT<=g_NAFVariable.KTABS;IT++)
2692	{
2693	/v0.96 M. GASTINEAU 18/12/98 : modification /
2694	/* TWIN[IT]=1.E0; //remplacee par:*/
2695	p_pardTWIN[IT]=1.E0;
2696	/v0.96 M. GASTINEAU 18/12/98 : fin modification /
2697	}
2698	}
2699	else if(g_NAFVariable.IW>=0)
2700	{
2701	CN = 1.E0;
2702	for(I = 1;I<=g_NAFVariable.IW;I++)
2703	{
2704	CN = I(2.E0/((double)(g_NAFVariable.IW+I))); /CN = CN2.D0I1.D0/(g_NAFVariable.IW+I)*/
2705	};
2706	for(IT=0;IT<=g_NAFVariable.KTABS;IT++)
2707	{
2708	T=IT*g_NAFVariable.XH-TM;
2709	/v0.96 M. GASTINEAU 18/12/98 : modification /
2710	/TWIN[IT]=CNpow((1.E0+cos(TPIST)),g_NAFVariable.IW);/
2711	/remplacee par:/
2712	p_pardTWIN[IT]=CNpow((1.E0+cos(TPIST)),g_NAFVariable.IW);
2713	/v0.96 M. GASTINEAU 18/12/98 : fin modification /
2714	}
2715	}
2716	else if(g_NAFVariable.IW==-1)
2717	{
2718	/v0.96 M. GASTINEAU 18/12/98 : modification /
2719	/* TWIN[0] =0.E0;
2720	TWIN[g_NAFVariable.KTABS] =0.E0;*/
2721	/remplacee par:/
2722	p_pardTWIN[0] =0.E0;
2723	p_pardTWIN[g_NAFVariable.KTABS] =0.E0;
2724	/v0.96 M. GASTINEAU 18/12/98 : fin modification /
2725	for(IT=1; IT<=g_NAFVariable.KTABS-1; IT++)
2726	{
2727	T=(IT*g_NAFVariable.XH-TM)/TM;
2728	/v0.96 M. GASTINEAU 18/12/98 : modification /
2729	/* TWIN[IT]= exp(-1.E0/(1.E0-TT))/CE;/
2730	/remplacee par:/
2731	p_pardTWIN[IT]= exp(-1.E0/(1.E0-T*T))/CE;
2732	/v0.96 M. GASTINEAU 18/12/98 : fin modification /
2733	}
2734	}
2735	if (g_NAFVariable.IPRT==1)
2736	{
2737	/!---------------- IMPRESSION TEMOIN/
2738	for(IT=0; IT<=g_NAFVariable.KTABS; (g_NAFVariable.KTABS>20?IT+=g_NAFVariable.KTABS/20:IT++))
2739	{
2740	/v0.96 M. GASTINEAU 18/12/98 : modification /
2741	/fprintf(g_NAFVariable.NFPRT,"%20.3E\n", TWIN[IT]1.E6);*/
2742	/remplacee par:/
2743	fprintf(g_NAFVariable.NFPRT,"%20.3E\n", p_pardTWIN[IT]*1.E6);
2744	/v0.96 M. GASTINEAU 18/12/98 : fin modification /
2745	}
2746	}
2747	}/* END SUBROUTINE INIWIN*/
2748
2749	/* SUBROUTINE ZARDYD(ZT,N,H,ZOM)*/
2750	BOOL naf_zardyd(t_complexe ZT, int N, double H, t_complexe ZOM)
2751	{
2752	/!**********************************************************************
2753	! CALCULE L'INTEGRALE D'UNE FONCTION TATULEE PAR LA METHODE DE HARDY
2754	! T(N) TABLEAU DOUBLE PRECISION DES VALEURS DE LA FONCTION
2755	! N = 6*K ENTIER
2756
2757	! H PAS ENTRE DEUX VALEURS
2758	! SOM VALEUR DE L'INTEGRALE SUR L'INTERVALLE [X1,XN]
2759	! LE PROGRAMME EST EN DOUBLE PRECISION
2760	! REVU LE 26/9/87 J. LASKAR
2761	!
2762	!************************************************************************
2763	IMPLICIT NONE
2764	! (ZT,N,H,ZOM)
2765	integer :: N
2766	complex (8) :: ZT(0:N),ZOM
2767	REAL (8) :: H
2768	!
2769	integer :: ITEST,K,I
2770	! */
2771	int ITEST, K, I;
2772	/* v0.96 M. GASTINEAU 13/01/99 : optimisation */
2773	#if NAF_USE_OPTIMIZE==0
2774	int INC;
2775	#else /NAF_USE_OPTIMIZE/
2776	double zomreel, zomimag;
2777	double zomreeltemp, zomimagtemp;
2778	double pdZT=(double)(ZT);
2779	#endif /NAF_USE_OPTIMIZE/
2780	/* v0.96 M. GASTINEAU 13/01/99 : fin optimisation */
2781
2782	ITEST=N%6; /ITEST=MOD(N,6);/
2783	if (ITEST!=0)
2784	{
2785	Myyerror("naf_zradyd - N N'EST PAS UN MULTIPLE DE 6\n");
2786	return FALSE;
2787	}
2788	K=N/6;
2789	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
2790	#if NAF_USE_OPTIMIZE==0
2791	/ZOM=41ZT(1)+216ZT(2)+27ZT(3)+272*ZT(4)
2792	$ +27ZT(5)+216ZT(6)+41ZT(N)/
2793	*ZOM = addcomplexe(
2794	addcomplexe(
2795	addcomplexe( muldoublcomplexe(41,ZT[0]), muldoublcomplexe(216,ZT[1]) )
2796	, addcomplexe( muldoublcomplexe(27,ZT[2]), muldoublcomplexe(272,ZT[3]) ) )
2797	, addcomplexe(
2798	addcomplexe( muldoublcomplexe(27,ZT[4]), muldoublcomplexe(216,ZT[5]) )
2799
2800	, muldoublcomplexe(41,ZT[N]) )
2801	);
2802	INC=0;
2803	for(I=1; I<=K-1; I++)
2804	{
2805	/* ZOM=ZOM+82ZT(6I+1)+216ZT(6I+2)+27ZT(6I+3)+272ZT(6I+4)
2806	$ +27ZT(6I+5)+216ZT(6I+6)*/
2807	INC +=6;
2808	*ZOM = addcomplexe(
2809	addcomplexe(
2810	addcomplexe( muldoublcomplexe(82,ZT[INC]), muldoublcomplexe(216,ZT[INC+1]) )
2811	, addcomplexe( muldoublcomplexe(27,ZT[INC+2]), muldoublcomplexe(272,ZT[INC+3]) ) )
2812	, addcomplexe(
2813	addcomplexe( muldoublcomplexe(27,ZT[INC+4]), muldoublcomplexe(216,ZT[INC+5]) )
2814	, *ZOM )
2815	);
2816	}
2817	ZOM= muldoublcomplexe(H6.E0/840.E0, ZOM); /ZOM =ZOMH6.D0/840.D0*/
2818	#else /remplacee par:/
2819	/* v0.96 M. GASTINEAU 13/01/99 : optimisation */
2820	#if 0
2821	/ZOM=41ZT(1)+216ZT(2)+27ZT(3)+272*ZT(4)
2822	$ +27ZT(5)+216ZT(6)+41ZT(N)/
2823	*ZOM = i_compl_add(
2824	i_compl_add(
2825	i_compl_add( i_compl_muldoubl(41,ZT[0]), i_compl_muldoubl(216,ZT[1]) )
2826	, i_compl_add( i_compl_muldoubl(27,ZT[2]), i_compl_muldoubl(272,ZT[3]) ) )
2827	, i_compl_add(
2828	i_compl_add( i_compl_muldoubl(27,ZT[4]), i_compl_muldoubl(216,ZT[5]) )
2829	, i_compl_muldoubl(41,ZT[N]) )
2830	);
2831	INC=0;
2832	for(I=1; I<=K-1; I++)
2833	{
2834	/* ZOM=ZOM+82ZT(6I+1)+216ZT(6I+2)+27ZT(6I+3)+272ZT(6I+4)
2835	$ +27ZT(6I+5)+216ZT(6I+6)*/
2836	INC +=6;
2837	*ZOM = i_compl_add(
2838	i_compl_add(
2839	i_compl_add( i_compl_muldoubl(82,ZT[INC]), i_compl_muldoubl(216,ZT[INC+1]) )
2840	, i_compl_add( i_compl_muldoubl(27,ZT[INC+2]), i_compl_muldoubl(272,ZT[INC+3]) ) )
2841	, i_compl_add(
2842	i_compl_add( i_compl_muldoubl(27,ZT[INC+4]), i_compl_muldoubl(216,ZT[INC+5]) )
2843	, *ZOM )
2844	);
2845	}
2846	ZOM= i_compl_muldoubl(H6.E0/840.E0, ZOM); /ZOM =ZOMH6.D0/840.D0*/
2847	#endif /0/ /remplacee par:/
2848	/ZOM=41ZT(1)+216ZT(2)+27ZT(3)+272*ZT(4)
2849	$ +27ZT(5)+216ZT(6)+41ZT(N)/
2850	zomreel=41 * ( pdZT++); zomimag=41 ( *pdZT++);
2851	zomreel+=216 * ( pdZT++); zomimag+=216 ( *pdZT++);
2852	zomreel+=27 * ( pdZT++); zomimag+=27 ( *pdZT++);
2853	zomreel+=272 * ( pdZT++); zomimag+=272 ( *pdZT++);
2854	zomreel+=27 * ( pdZT++); zomimag+=27 ( *pdZT++);
2855	zomreel+=216 * ( pdZT++); zomimag+=216 ( *pdZT++);
2856	zomreel+=41 * ( ZT[N].reel); zomimag+=41 * ( ZT[N].imag);
2857
2858	for(I=1; I<=K-1; I++)
2859	{
2860	/* ZOM=ZOM+82ZT(6I+1)+216ZT(6I+2)+27ZT(6I+3)+272ZT(6I+4)
2861	$ +27ZT(6I+5)+216ZT(6I+6)*/
2862	zomreeltemp=82 * ( pdZT++); zomimagtemp=82 ( *pdZT++);
2863	zomreeltemp+=216 * ( pdZT++); zomimagtemp+=216 ( *pdZT++);
2864	zomreeltemp+=27 * ( pdZT++); zomimagtemp+=27 ( *pdZT++);
2865	zomreeltemp+=272 * ( pdZT++); zomimagtemp+=272 ( *pdZT++);
2866	zomreeltemp+=27 * ( pdZT++); zomimagtemp+=27 ( *pdZT++);
2867	zomreeltemp+=216 * ( pdZT++); zomimagtemp+=216 ( *pdZT++);
2868	zomreel += zomreeltemp; zomimag += zomimagtemp;
2869	}
2870	ZOM->reel = H6.E0/840.E0zomreel;
2871	ZOM->imag = H6.E0/840.E0zomimag;
2872	/* v0.96 M. GASTINEAU 13/01/99 : fin d'optimisation */
2873	#endif /NAF_USE_OPTIMIZE==0/
2874	/* v0.96 M. GASTINEAU 01/12/98 : fin optimisation */
2875	return TRUE;
2876	}/* END SUBROUTINE ZARDYD*/
2877
2878
2879	/* SUBROUTINE PROSCAA(F1,F2,ZP)*/
2880	BOOL naf_proscaa(double F1, double F2, t_complexe *ZP)
2881	{
2882	/*!-----------------------------------------------------------------------
2883	! PROSCAA CALCULE LE PRODUIT SCALAIRE
2884	! < EXP(IF1T), EXP(IF2T)>
2885	! SUR L'INTERVALLE [0:g_NAFVariable.KTABS] T=g_NAFVariable.T0+g_NAFVariable.XH*IT
2886	! CALCUL NUMERIQUE
2887	! ZP PRODUIT SCALAIRE COMPLEXE
2888	! F1,F2 FREQUENCES EN "/AN
2889	! REVU LE 26/9/87 J. LASKAR
2890	!-----------------------------------------------------------------------
2891
2892	IMPLICIT NONE
2893	! (g_NAFVariable.KTABS,F1,F2,ZP)
2894	REAL (8) :: F1,F2
2895	complex (8) :: ZP
2896	!
2897	integer :: LTF
2898	real *8 :: OM,ANG0,ANGI,H
2899	complex (8) :: ZI,ZAC,ZINC,ZEX
2900	complex (8), dimension (:),allocatable :: ZTF
2901
2902	!*/
2903	int LTF;
2904	double OM,ANG0,ANGI,H;
2905	t_complexe ZI,ZAC,ZINC,ZEX;
2906	t_complexe *ZTF=NULL;
2907
2908	SYSCHECKMALLOCSIZE(ZTF, t_complexe, g_NAFVariable.KTABS+1); /allocate(ZTF(0:g_NAFVariable.KTABS),stat = nerror)/
2909	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
2910	#if NAF_USE_OPTIMIZE==0
2911	ZI=cmplx(0.E0,1.E0);
2912	/!----------! FREQUENCES EN UNITE D'ANGLE PAR UNITE DE TEMPS/
2913	OM = (F1-F2)/g_NAFVariable.UNIANG;
2914	LTF=g_NAFVariable.KTABS;
2915	ANG0=OM*g_NAFVariable.T0;
2916	ANGI=OM*g_NAFVariable.XH;
2917	ZAC= expcomplexe(muldoublcomplexe(-ANG0, ZI)); /ZAC = EXP (-ZIANG0)*/
2918	ZINC=expcomplexe(muldoublcomplexe(-ANGI, ZI)); /ZINC= EXP (-ZIANGI)*/
2919	ZEX= divcomplexe(ZAC, ZINC); /ZEX = ZAC/ZINC/
2920	naf_ztpow2a(g_NAFVariable.KTABS,64,ZTF,TWIN,ZINC,ZEX);/CALL ZTPOW2A(g_NAFVariable.KTABS+1,64,ZTF,TWIN,ZINC,ZEX)/
2921	#else/remplacee par:/
2922	i_compl_cmplx(&ZI,0.E0,1.E0);
2923	/!----------! FREQUENCES EN UNITE D'ANGLE PAR UNITE DE TEMPS/
2924	OM = (F1-F2)/g_NAFVariable.UNIANG;
2925	LTF=g_NAFVariable.KTABS;
2926	ANG0=OM*g_NAFVariable.T0;
2927	ANGI=OM*g_NAFVariable.XH;
2928	ZAC= i_compl_exp(i_compl_muldoubl(-ANG0, ZI)); /ZAC = EXP (-ZIANG0)*/
2929	ZINC=i_compl_exp(i_compl_muldoubl(-ANGI, ZI)); /ZINC= EXP (-ZIANGI)*/
2930	ZEX= i_compl_div(ZAC, ZINC); /ZEX = ZAC/ZINC/
2931	naf_ztpow2a(g_NAFVariable.KTABS,64,ZTF,TWIN,ZINC,ZEX);/CALL ZTPOW2A(g_NAFVariable.KTABS+1,64,ZTF,TWIN,ZINC,ZEX)/
2932	#endif /NAF_USE_OPTIMIZE==0/
2933	/* v0.96 M. GASTINEAU 01/12/98 : fin optimisation */
2934
2935	/!------------------ TAILLE DU PAS/
2936	H=1.E0/LTF;
2937	/CALL ZARDYD(ZTF,LTF+1,H,ZP)/
2938	if (naf_zardyd(ZTF,LTF,H,ZP)==FALSE)
2939	{
2940	SYSFREE(ZTF);
2941	return FALSE;
2942	}
2943	SYSFREE(ZTF);
2944	return TRUE;
2945	}/* END SUBROUTINE PROSCAA*/
2946
2947	/* SUBROUTINE ZTPOW2A (N,N1,ZTF,TW,ZA,ZAST)*/
2948	void naf_ztpow2a(int N, int N1, t_complexe ZTF, double TW, t_complexe ZA, t_complexe ZAST)
2949	{
2950	/*!-----------------------------------------------------------------------
2951	! ZTPOW CALCULE ZTF(I) = TW(I)ZASTZA**I EN VECTORIEL
2952	! ZT(N1) ZONE DE TRAVAIL
2953	!-----------------------------------------------------------------------
2954	IMPLICIT NONE
2955	! (N,N1,ZT,ZTF,TW,ZA,ZAST)
2956	integer :: N,N1
2957	real *8 :: TW(0:N)
2958	complex (8) :: ZTF(0:N),ZA,ZAST
2959	!
2960	integer :: I,INC,IT, NX,NT
2961	complex (8) :: ZT1,ZINC
2962	complex (8), dimension (:),allocatable :: ZT
2963	! */
2964	int I,INC,IT, NX,NT;
2965	t_complexe ZT1,ZINC;
2966	t_complexe *ZT=NULL;
2967
2968	SYSCHECKMALLOCSIZE(ZT, t_complexe, N1); /allocate(ZT(0:N1-1),stat = nerror)/
2969	if (N<N1-1)
2970	{
2971	fprintf(stdout,"DANS ZTPOW, N = %d\n",N);
2972	return;
2973	}
2974	/!----------! /
2975	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
2976	#if NAF_USE_OPTIMIZE==0
2977	ZT[0] = mulcomplexe(ZAST,ZA);
2978	for(I = 1; I<N1; I++)
2979	{
2980	ZT[I] = mulcomplexe(ZT[I-1],ZA);
2981	}
2982	for(I = 0;I<N1; I++)
2983	{
2984	ZTF[I] = muldoublcomplexe(TW[I],ZT[I]);
2985	}
2986	ZT1 = divcomplexe(ZT[N1-1],ZAST); /ZT1 = ZT(N1-1)/ZAST/
2987	ZINC= cmplx(1E0,0E0);
2988	INC =0;
2989	NT = (N+1)/N1 ;
2990	for(IT = 2;IT<= NT; IT++)
2991	{
2992	ZINC=mulcomplexe(ZINC, ZT1); /ZINC = ZINCZT1*/
2993	INC += N1; /INC = INC + N1/
2994	for(I = 0; I<N1; I++)
2995	{
2996	ZTF[INC +I] = muldoublcomplexe(TW[INC+I],mulcomplexe(ZT[I],ZINC));/ZTF(INC +I) = TW(INC+I)ZT(I)ZINC/
2997
2998	}
2999	}
3000	ZINC=mulcomplexe(ZINC, ZT1); /ZINC = ZINCZT1*/
3001	INC += N1; /INC = INC + N1/
3002	NX = N+1-NT*N1;
3003	for(I = 0;I< NX; I++)
3004	{
3005	ZTF[INC +I] = muldoublcomplexe(TW[INC+I],mulcomplexe(ZT[I],ZINC));/ZTF(INC +I) = TW(INC+I)ZT(I)ZINC/
3006	}
3007	#else /remplacee par:/
3008	ZT[0] = i_compl_mul(ZAST,ZA);
3009	for(I = 1; I<N1; I++)
3010	{
3011	ZT[I] = i_compl_mul(ZT[I-1],ZA);
3012	}
3013	for(I = 0;I<N1; I++)
3014	{
3015	ZTF[I] = i_compl_muldoubl(TW[I],ZT[I]);
3016	}
3017	ZT1 = i_compl_div(ZT[N1-1],ZAST); /ZT1 = ZT(N1-1)/ZAST/
3018	i_compl_cmplx(&ZINC,1E0,0E0);
3019	INC =0;
3020	NT = (N+1)/N1 ;
3021	for(IT = 2;IT<= NT; IT++)
3022	{
3023	i_compl_pmul(&ZINC,&ZT1); /ZINC = ZINCZT1*/
3024	INC += N1; /INC = INC + N1/
3025	for(I = 0; I<N1; I++)
3026	{
3027	ZTF[INC +I] = i_compl_muldoubl(TW[INC+I],i_compl_mul(ZT[I],ZINC));/ZTF(INC +I) = TW(INC+I)ZT(I)ZINC/
3028
3029	}
3030	}
3031	i_compl_pmul(&ZINC,&ZT1); /ZINC = ZINCZT1*/
3032	INC += N1; /INC = INC + N1/
3033	NX = N+1-NT*N1;
3034	for(I = 0;I< NX; I++)
3035	{
3036	ZTF[INC +I] = i_compl_muldoubl(TW[INC+I],i_compl_mul(ZT[I],ZINC));/ZTF(INC +I) = TW(INC+I)ZT(I)ZINC/
3037	}
3038	#endif /NAF_USE_OPTIMIZE==0/
3039	/* v0.96 M. GASTINEAU 01/12/98 : fin optimisation */
3040	SYSFREE(ZT);
3041	}/* END SUBROUTINE ZTPOW2A*/
3042
3043	/* SUBROUTINE CORRECTION(FREQ)*/
3044	void naf_correction(double *FREQ)
3045	{
3046	/*!------------------------------------------------------------------
3047	! RETOURNE LA PREMIERE FREQUENCE (FREQ) CORRIGE PAR LES SUIVANTES
3048	!
3049	! F. Joutel 26/5/98
3050	! Ref : J. Laskar : Introduction to frequency map analysis
3051	!------------------------------------------------------------------
3052
3053	IMPLICIT NONE
3054	! (FREQ)
3055	real *8 :: FREQ
3056	!
3057	integer :: I
3058	real *8 :: TL,TM2,PICARRE,TM,FAC,A,COR,OMEGA,DELTA
3059	complex (8) :: ZI,ZALPHA
3060	!
3061	! PI = ATAN2(1.D0,1.D0)4.D0/
3062	int I;
3063
3064	double TL,TM2,PICARRE,TM,FACTEUR,A,COR,OMEGA,DELTA;
3065	t_complexe ZI,ZALPHA;
3066
3067	TL=g_NAFVariable.KTABSg_NAFVariable.XH0.5E0;
3068	TM2=1.E0/(TL*TL);
3069	PICARRE=pi*pi;
3070	/! g_NAFVariable.IW >= 1/
3071	TM=g_NAFVariable.T0+TL;
3072	i_compl_cmplx(&ZI,0.E0,1.E0);
3073	if (g_NAFVariable.IW>=0)
3074	{
3075	FACTEUR = -TM2;
3076	A=-1.E0/3.E0;
3077	for(I=1;I<=g_NAFVariable.IW;I++)
3078	{
3079	FACTEUR = -PICARRETM2(II); /FACTEUR=-FACTEURPICARRETM2(I*2)/
3080	A+=2.E0/(PICARREII); /A=A+2.D0/(PICARREI*2) /
3081	}
3082	FACTEUR /=A; /FACTEUR = FACTEUR/A/
3083	COR=0.E0;
3084	for(I=2;I<=g_NAFVariable.NFS;I++)
3085	{
3086	OMEGA=(g_NAFVariable.TFS[I]-g_NAFVariable.TFS[1])/g_NAFVariable.UNIANG;
3087	/* v0.96 M. GASTINEAU 01/12/98 : utilisation des fonctions complexes inlines et optimisation */
3088	#if NAF_USE_OPTIMIZE==0
3089	ZALPHA = mulcomplexe(divcomplexe(g_NAFVariable.ZAMP[I],g_NAFVariable.ZAMP[1]), expcomplexe(muldoublcomplexe(OMEGATM,ZI))); /ZALPHA=g_NAFVariable.ZAMP[I]/g_NAFVariable.ZAMP[1]EXP(ZIOMEGATM)/
3090	#else /remplacee par:/
3091	ZALPHA = i_compl_mul(i_compl_div(g_NAFVariable.ZAMP[I],g_NAFVariable.ZAMP[1]),
3092	i_compl_exp(i_compl_muldoubl(OMEGA*TM,ZI)));
3093	/ZALPHA=g_NAFVariable.ZAMP[I]/g_NAFVariable.ZAMP[1]EXP(ZIOMEGATM)*/
3094	#endif /NAF_USE_OPTIMIZE==0/
3095	/* v0.96 M. GASTINEAU 01/12/98 : fin optimisation */
3096	DELTA=FACTEURZALPHA.reel/(pow(OMEGA,(2g_NAFVariable.IW+1)))cos(OMEGATL); /DELTA=FACTEURDREAL(ZALPHA)/OMEGA*(2g_NAFVariable.IW+1)COS(OMEGATL)*/
3097	COR += DELTA; /COR=COR+DELTA/
3098	if (g_NAFVariable.IPRT>=2)
3099	{
3100	fprintf(g_NAFVariable.NFPRT,"OMEGA=%g, DELTA=%g ,COR=%g\n",OMEGA,DELTA,COR);
3101	}
3102	}
3103	COR = g_NAFVariable.UNIANG; /COR=CORg_NAFVariable.UNIANG/
3104	*FREQ=g_NAFVariable.TFS[1]+COR;
3105	if (g_NAFVariable.IPRT>=1)
3106	{
3107	fprintf(g_NAFVariable.NFPRT,"CORRECTION DE LA 1ERE FREQUENCE\n");
3108	fprintf(g_NAFVariable.NFPRT, "FREQ. DE DEPART=%20.15E, CORRECTION=%20.15E\n",g_NAFVariable.TFS[1],COR);
3109	fprintf(g_NAFVariable.NFPRT, "FREQUENCE CORRIGEE:%20.15E\n", *FREQ);
3110	}
3111	}
3112	else
3113	{
3114	if (g_NAFVariable.IPRT>=0)
3115	{
3116	fprintf(g_NAFVariable.NFPRT, "PAS DE CALCUL\n");
3117	}
3118	}
3119	}/* END SUBROUTINE CORRECTION*/
3120
3121	/* end module naff
3122	*/
3123	/-----------------------------------------------------------/
3124	/* cree et retourneun element de liste de fenetre pour naf. */
3125	/* Les champs sont remplis avec les arguments de la fonction.*/
3126	/* L'element retourne doit etre libere avec */
3127	/* delete_list_fenetre_naf. */
3128	/-----------------------------------------------------------/
3129	/* v0.96 M. GASTINEAU 18/12/98 : ajout */
3130	t_list_fenetre_naf *cree_list_fenetre_naf(const double p_dFreqMin, const double p_dFreqMax,
3131	const int p_iNbTerm)
3132	{
3133	t_list_fenetre_naf *pListRes;
3134	#ifdef DBUG
3135	if (niveau>=4)
3136	{
3137	printf("cree_list_fenetre_naf(%g, %g, %d) - entree \n",
3138	p_dFreqMin, p_dFreqMax, p_iNbTerm);
3139	}
3140	#endif /DBUG/
3141	SYSCHECKMALLOC(pListRes, t_list_fenetre_naf);
3142	pListRes->suivant = NULL;
3143	pListRes->dFreqMin = p_dFreqMin;
3144	pListRes->dFreqMax = p_dFreqMax;
3145	pListRes->iNbTerme = p_iNbTerm;
3146	#ifdef DBUG
3147	if (niveau>=4)
3148	{
3149	printf("cree_list_fenetre_naf() - sortie et retourne 0x%p \n",pListRes);
3150	}
3151	#endif /DBUG/
3152	return pListRes;
3153	}
3154
3155	/-----------------------------------------------------------/
3156	/Detruit une liste cree avec cree_list_fenetre_naf. /
3157	/* En sortie, p_pListFenNaf pointe vers une adresse invalide.*/
3158	/-----------------------------------------------------------/
3159	/* v0.96 M. GASTINEAU 18/12/98 : ajout */
3160	void delete_list_fenetre_naf(t_list_fenetre_naf *p_pListFenNaf)
3161	{
3162	#ifdef DBUG
3163	if (niveau>=4)
3164	printf("delete_list_fenetre_naf(0x%p) - entree \n",p_pListFenNaf);
3165	#endif /DBUG/
3166	while(p_pListFenNaf!=NULL)
3167	{
3168	t_list_fenetre_naf *next=p_pListFenNaf->suivant;
3169	SYSFREE(p_pListFenNaf);
3170	p_pListFenNaf=next;
3171	}
3172
3173	#ifdef DBUG
3174	if (niveau>=4)
3175	printf("delete_list_fenetre_naf() - sortie \n");
3176	#endif /DBUG/
3177	}
3178
3179	/-----------------------------------------------------------/
3180	/* Concatene 2 listes de fenetre de naf. */
3181	/* En sortie, p_pListFenHead et p_pListFenEnd ne doivent pas */
3182	/* etre liberee. Seul la liste retournee doit etre liberee. */
3183	/* Les deux parametres peuvent etre NULL. */
3184	/-----------------------------------------------------------/
3185	/* v0.96 M. GASTINEAU 18/12/98 : ajout */
3186	t_list_fenetre_naf concat_list_fenetre_naf(t_list_fenetre_naf p_pListFenHead,
3187	t_list_fenetre_naf *p_pListFenEnd)
3188	{
3189	t_list_fenetre_naf *pListTemp;
3190	#ifdef DBUG
3191	if (niveau>=4)
3192	printf("concat_list_fenetre_naf(0x%p, 0x%p) - entree \n",p_pListFenHead,p_pListFenEnd);
3193	#endif /DBUG/
3194	if (p_pListFenHead==NULL)
3195	{
3196	return p_pListFenEnd;
3197	}
3198	else if (p_pListFenEnd==NULL)
3199	{
3200	return p_pListFenHead;
3201	}
3202	else
3203	{ /* p_pListFenHead is not NULL */
3204	for(pListTemp=p_pListFenHead;
3205	pListTemp->suivant!=NULL;
3206	pListTemp = pListTemp->suivant)
3207	{ /do nothing /}
3208	/concatene /
3209	pListTemp->suivant = p_pListFenEnd;
3210	return p_pListFenHead;
3211	}
3212	}

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: TRACY3/trunk/tracy/tracy/src/modnaff.cc @ 3

Download in other formats: