Context Navigation

source: MML/trunk/machine/SOLEIL/common/naff/nafflib/modnaff.c @ 4

Last change on this file since 4 was 4, checked in by zhangj, 10 years ago
Initial import--MML version from SOLEIL@2013
File size: 120.8 KB

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

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

Download in other formats: