?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fix saving lists of arrays with recent versions of numpy
[qpms.git] / qpms / gaunt.c
blob33c4da7f78296a6f7a5df25b6b5c75538b50c8d7
1 #include "gaunt.h"
2 #ifdef USE_FORTRAN_GAUNT_XU
3 void __vec_trans_MOD_gaunt_xu(const double *m, const double *n, const double *mu, const double *nu, const int *qmax, double *v_aq, int *err);
4
5 void gaunt_xu(int m, int n, int mu, int nu, int qmax, double *v_aq, int *err) {
6 double mf = m, nf=n, muf=mu,nuf=nu;
7 __vec_trans_MOD_gaunt_xu(&mf,&nf,&muf,&nuf,&qmax,v_aq,err);
8 }
9
10 #else //!USE_FORTRAN_GAUNT_XU
11 //#include <stdlib.h>
12 #include <math.h>
13 #include <stdio.h>
14 #include "assert_cython_workaround.h"
15
16 #define ZERO_THRESHOLD 1.e-8
17 #define BF_PREC 1.e-10
18 // "Besides, the determined Real Programmer can write FORTRAN programs in any language."
19 // -- Ed Post, Real Programmers Don't Use Pascal, 1982.
20
21 // logarithm of factorial (from basicsubs.f90)
22 double lnf (double z) {
23 // expansion parameters
24 static const double v_c0[] = {
25 0.16427423239836267e5, -0.48589401600331902e5, 0.55557391003815523e5, -0.30964901015912058e5,
26 0.87287202992571788e4, -0.11714474574532352e4, 0.63103078123601037e2, -0.93060589791758878e0,
27 0.13919002438227877e-2,-0.45006835613027859e-8, 0.13069587914063262e-9
28 };
29 static const double cp = 2.5066282746310005;
30 double a = 1.;
31 double b = z + 10.5;
32 b = (z + 0.5) * log(b) - b;
33 for (int i = 0; i < (sizeof(v_c0) / sizeof(double)); i++) {
34 z += 1.;
35 a += v_c0[i] / z;
36 }
37
38 return b+log(cp*a);
39 }
40
41 // logarithm of Pochhammer function (from basicsubs.f90)
42 // FIXME replace with standard C99 lgamma functions
43 double lpoch(double x, double n) {
44 if(fabs(n) < 1e-5) // ???
45 return 1.;
46 double sum = x+n;
47 return lnf(sum-1.) - lnf(x-1.);
48 }
49
50 // pochhammer function: substitute for fortran DPOCH
51 // gamma(a+x) / gamma(a)
52 double poch(double a, double x) { // FIXME replace with standard C99 lgamma functions
53 return exp(lpoch(a,x));
54 }
55
56 double f_a0 (int m, int n, int mu, int nu) {
57 double logw = lnf(n+nu-m-mu) - lnf(n-m) - lnf(nu-mu);
58 double logp = lpoch(n+1, n) + lpoch(nu+1, nu) - lpoch(n+nu+1, n+nu);
59 return exp(logw+logp);
60 }
61
62 // coefficient Ap(m,n,mu,nu,p) (from vec_trans.f90)
63 int f_Ap(int m, int n, int mu, int nu, int p) {
64 return p*(p-1)*(m-mu)-(m+mu)*(n-nu)*(n+nu+1);
65 }
66
67 // coefficient a(m,n,mu,nu,1) normalized by the backward recursion (from vec_trans.f90)
68
69 double f_a1norm(int m, int n, int mu, int nu) {
70 int n4 = n + nu - m - mu;
71 return ((2.*n + 2.*nu - 3.) / 2.) * (1. - ((2.*n + 2.*nu - 1.) / (n4 * (n4-1.)))
72 * ((m - n) * (m - n + 1.) / (2.*n - 1.) + (mu-nu) * (mu-nu+1.)/(2.*nu-1.)));
73 }
74
75 // coefficient a(m,n,mu,nu,2) normalized by the backward recursion (From vec_trans.f90)
76
77 double f_a2norm(int m, int n, int mu, int nu) {
78 double n4 = n + nu - m - mu;
79 double n2nu = 2*n + 2*nu;
80 double mn = m - n;
81 double munu = mu - nu;
82
83 return ((n2nu-1.)*(n2nu-7.)/4.) * \
84 ( ((n2nu-3.)/(n4*(n4-1.))) * \
85 ( ((n2nu-5.)/(2.*(n4-2.)*(n4-3.))) * \
86 ( mn*(mn+1.)*(mn+2.)*(mn+3.)/((2.*n-1.)*(2.*n-3.)) + \
87 2.*mn*(mn+1.)*munu*(munu+1.)/((2.*n-1.)*(2.*nu-1.)) + \
88 munu*(munu+1.)*(munu+2.)*(munu+3.)/((2.*nu-1.)*(2.*nu-3.)) \
89 ) - mn*(mn+1.)/(2.*n-1.) - munu*(munu+1.)/(2.*nu-1.) ) +0.5);
90 }
91
92 // just for convenience – square of ints
93 static inline int isq(int x) {return x * x;}
94 static inline double fsq(float x) {return x * x;}
95
96 double f_alpha(int n, int nu, int p) {
97 return (isq(p) - isq(n-nu))*(isq(p)-isq(n+nu+1)) / (double)(4*isq(p)-1);
98 }
99
100 static inline int min1pow(int pow) { return (pow % 2) ? -1 : 1; }
101
102 // starting value of coefficient a(m,n,mu,nu,qmax) for the forward recursion
103 double f_aqmax(int m, int n, int mu, int nu, int qmax) {
104 int pmin = n + nu - 2*qmax;
105 // TODO? zde měl int a double varianty téhož – má to nějaký smysl???
106 if (pmin == n-nu) {
107 double logw = lnf(n+m) + lnf(2*pmin+1) - lnf(nu-mu) - lnf(n-nu) - lnf(pmin+m+mu);
108 double logp = lpoch(nu+1, nu) - lpoch(n+1, n+1);
109 return min1pow(mu)*exp(logw+logp);
110 } else if (pmin == nu-n) {
111 double logw = lnf(nu+mu) + lnf(2*pmin+1) - lnf(n-m) - lnf(nu-n) - lnf(pmin+m+mu);
112 double logp = lpoch(n+1, n) - lpoch(nu+1, nu+1); // ??? nešel by druhý lpoch nahradit něčím rozumnějším?
113 return min1pow(m)*exp(logw+logp);
114 } else if (pmin == m+mu) {
115 double logw = lpoch(qmax+1,qmax)+lnf(n+nu-qmax)+lnf(n+m)+lnf(nu+mu) \
116 -lnf(n-qmax)-lnf(nu-qmax)-lnf(n-m)-lnf(nu-mu)-lnf(n+nu+pmin+1);
117 return min1pow(n+m-qmax)*(2*pmin+1)*exp(logw);
118 } else if (pmin == -m-mu) {
119 double logw = lpoch(qmax+1,qmax)+lnf(n+nu-qmax)+lnf(pmin-m-mu) \
120 -lnf(n-qmax)-lnf(nu-qmax)-lnf(n+nu+pmin+1);
121 return min1pow(nu+mu-qmax)*(2*pmin+1)*exp(logw);
122 } else if (pmin==m+mu+1) {
123 int Apmin = f_Ap(m,n,mu,nu,pmin);
124 double logw = lpoch(qmax+1,qmax)+lnf(n+nu-qmax)+lnf(n+m)+lnf(nu+mu) \
125 -lnf(n+nu+pmin+1)-lnf(n-qmax)-lnf(nu-qmax)-lnf(n-m)-lnf(nu-mu);
126 return min1pow(n+m-qmax)*Apmin*(2*pmin+1)*exp(logw)/(double)(pmin-1);
127 } else if (pmin==-m-mu+1) {
128 int Apmin=f_Ap(m,n,mu,nu,pmin);
129 double logw=lpoch(qmax+1,qmax)+lnf(n+nu-qmax)+lnf(pmin-m-mu) \
130 -lnf(n+nu+pmin+1)-lnf(n-qmax)-lnf(nu-qmax);
131 return min1pow(nu+mu-qmax)*Apmin*(2*pmin+1)*exp(logw) / (double)(pmin-1);
132 }
133 assert(0);
134 }
135
136 // starting value of coefficient a(m,n,mu,nu,qmax-1) for the forward recursion
137 double f_aqmax_1(int m, int n, int mu, int nu, int qmax) {
138 int pmin=n+nu-2*qmax;
139 //pmin_i=INT(pmin,lo)
140 //qmaxr=REAL(qmax,dbl)
141 int Apmin2=f_Ap(m,n,mu,nu,pmin+2);
142
143 if (pmin==(m+mu+1)) { // pmin_if
144 double logw=lpoch(qmax+1,qmax)+lnf(n+nu-qmax)+lnf(n+m)+lnf(nu+mu)
145 -lnf(n+nu+pmin)-lnf(n-qmax+1)-lnf(nu-qmax+1)-lnf(n-m)-lnf(nu-mu);
146 double f1=min1pow(m+n-qmax)*Apmin2*(2.*pmin+3.)*(2.*pmin+5.) /(pmin*(n+nu+pmin+3.));
147 double f2=(n-qmax)*(nu-qmax)*(2.*qmax+1.)/((pmin+m+mu+1.)*(pmin+m+mu+2.)*(2*qmax-1.));
148 return f1*f2*exp(logw);
149 } else if (pmin==(-m-mu+1)) {
150 double logw=lpoch(qmax+1.,qmax+1.)+lnf(n+nu-qmax)+lnf(pmin-m-mu)
151 -lnf(n+nu+pmin)-lnf(n-qmax+1.)-lnf(nu-qmax+1.);
152 double f1=min1pow(nu+mu-qmax)*Apmin2*(2.*pmin+3.)*(2.*pmin+5.)
153 /(6.*pmin*(n+nu+pmin+3.));
154 double f2=(n-qmax)*(nu-qmax)/(2.*qmax-1.);
155 return f1*f2*exp(logw);
156 } // END IF pmin_if
157 assert(0);
158 }//END FUNCTION f_aqmax_1
159
160
161 // coeff a(m,n,mu,nu,2) normalizzato per la backward recursion calcolato questa volta per ricorsione
162 // (from vec_trans.f90)
163 double f_a2normr(int m, int n, int mu, int nu, double a1norm) {
164 int p = n + nu - 4;
165 int p1 = p - m - mu;
166 int p2 = p + m + mu;
167 int Ap4 = f_Ap(m,n,mu,nu,p+4);
168 int Ap6 = f_Ap(m,n,mu,nu,p+6);
169
170 double alphap1=f_alpha(n,nu,p+1);
171 double alphap2=f_alpha(n,nu,p+2);
172
173 if (!Ap4) {
174 if(!Ap6) {
175
176 double c0=(p+2.)*(p1+1.)*alphap1;
177 double c1=(p+1.)*(p2+2.)*alphap2;
178
179 return (c1/c0)*a1norm;
180 } else /* Ap6 != 0 */ {
181 double Ap2=f_Ap(m,n,mu,nu,p+2);
182 double Ap3=f_Ap(m,n,mu,nu,p+3);
183 double Ap5=f_Ap(m,n,mu,nu,p+5);
184 double alphap5=f_alpha(n,nu,p+5);
185
186 double c0=(p+2.)*(p+3.)*(p+5.)*(p1+1.)*(p1+2.)*(p1+4.)*Ap6*alphap1;
187 double c1=(p+5.)*(p1+4.)*Ap6*(Ap2*Ap3+(p+1.)*(p+3.)*(p1+2.)*(p2+2.)*alphap2);
188 double c2=(p+2.)*(p2+3.)*Ap2*(Ap5*Ap6+(p+4.)*(p+6.)*(p1+5.)*(p2+5.)*alphap5);
189
190 return (c1/c0)*a1norm+(c2/c0);
191 }
192 } else /* Ap4 != 0 */ {
193 double Ap2=f_Ap(m,n,mu,nu,p+2);
194 double Ap3=f_Ap(m,n,mu,nu,p+3);
195 double alphap3=f_alpha(n,nu,p+3);
196 double alphap4=f_alpha(n,nu,p+4);
197
198 double c0=(p+2.)*(p+3.)*(p1+1.)*(p1+2.)*Ap4*alphap1;
199 double c1=Ap2*Ap3*Ap4+(p+1.)*(p+3.)*(p1+2.)*(p2+2.)*Ap4*alphap2 \
200 + (p+2.)*(p+4.)*(p1+3.)*(p2+3.)*Ap2*alphap3;
201 double c2=-(p+2.)*(p+3.)*(p2+3.)*(p2+4.)*Ap2*alphap4;
202
203 return (c1/c0)*a1norm+(c2/c0);
204 }
205 }
206
207 #define MAX(x,y) (((x) > (y)) ? (x) : (y))
208
209 int j3minimum(int j1, int j2, int m1, int m2) {
210 return MAX(abs(j1-j2), abs(m1+m2));
211 }
212
213 int nw(int j1, int j2, int m1, int m2) {
214 return j1+j2+1-MAX(abs(j1-j2),abs(m1+m2));
215 }
216
217 // per calcolare il primo termine della backward recursion; vec_trans.f90
218 double wdown0(int j1, int j2, int m1, int m2) {
219 double logw=.5*(lnf(2.*j1)+lnf(2.*j2)+lnf(j1+j2-m1-m2)+lnf(j1+j2+m1+m2) - \
220 lnf(2.*j1+2.*j2+1.)-lnf(j1-m1)-lnf(j1+m1)-lnf(j2-m2)-lnf(j2+m2));
221 return min1pow(j1+j2+m1+m2)*exp(logw);
222 }
223
224 // per calcolar il primo termine della upward recursion; vec_trans.f90
225 double wup0(int j1, int j2, int m1, int m2) {
226 double logw;
227 if ( ((j1-j2)>=0) && ( abs(j1-j2)>=abs(m1+m2)) ) {
228
229 logw=0.5 * ( lnf(j1-m1)+lnf(j1+m1)+lnf(2.*j1-2.*j2)+lnf(2.*j2) - \
230 lnf(j2-m2)-lnf(j2+m2)-lnf(j1-j2-m1-m2)-lnf(j1-j2+m1+m2)-lnf(2.*j1+1.) );
231
232 return min1pow(j1+m1)*exp(logw);
233 } else if ( ((j1-j2)<0) && ( abs(j1-j2)>=abs(m1+m2)) ) {
234
235 logw=0.5 * ( lnf(j2-m2)+lnf(j2+m2)+lnf(2.*j2-2.*j1)+lnf(2.*j1) - \
236 lnf(j1-m1)-lnf(j1+m1)-lnf(j2-j1-m1-m2)-lnf(j2-j1+m1+m2)-lnf(2.*j2+1.) );
237
238 return min1pow(j2+m2)*exp(logw);
239 } else if ( ((m1+m2)>0) && (abs(j1-j2)<abs(m1+m2)) ) {
240
241 logw=0.5 * ( lnf(j1+m1)+lnf(j2+m2)+lnf(j1+j2-m1-m2)+lnf(2.*m1+2.*m2) - \
242 lnf(j1-m1)-lnf(j2-m2)-lnf(j1-j2+m1+m2)-lnf(j2-j1+m1+m2)-lnf(j1+j2+m1+m2+1.) );
243
244 return min1pow(j2+m2)*exp(logw);
245 } else if ( ((m1+m2)<0) && (abs(j1-j2)<abs(m1+m2)) ) {
246
247 logw=0.5 * ( lnf(j1-m1)+lnf(j2-m2)+lnf(j1+j2+m1+m2)+lnf(-2.*m1-2.*m2) - \
248 lnf(j1+m1)-lnf(j2+m2)-lnf(j1-j2-m1-m2)-lnf(j2-j1-m1-m2)-lnf(j1+j2-m1-m2+1.) );
249
250 return min1pow(j1+m1)*exp(logw);
251 }
252 assert(0);
253 }
254
255 // coefficiente per il calcolo ricorsivo
256 double cr(double j1, double j2, double j3, double m1, double m2, double m3) {
257 return (fsq(j3)-fsq(j1-j2))*( fsq(j1+j2+1)-fsq(j3) )*( fsq(j3)-fsq(m3) );
258 }
259
260
261 // coefficiente per il calcolo ricorsivo
262 // (nebo raději všude doubly)?
263 double dr(double j1, double j2, double j3, double m1, double m2, double m3) {
264 return -(2*j3+1)*( j1*(j1+1)*m3-j2*(j2+1)*m3-j3*(j3+1)*(m2-m1) );
265 }
266
267
268 //******************************************************************************
269 //7) subroutine Wigner3jm: calcolo il vettore di simboli 3jm
270 //******************************************************************************
271
272 void wigner3jm(int j1, int j2, int m1, int m2, int j3min, int j3max, double *v_w3jm){
273 // in the original code, the dimension of v_w3jm is (j3min:j3max).
274 // In C, this means it has length j3max-j3min+1, and we must
275 // always deduct j3min from the indices
276
277 // Inizializzo gli indici per la downward recursion
278 int j3 = j3max; // we don't use separate j3int as gevero does.
279
280 // In questo if separo i casi in cui ho un vettore di lunghezza uno da quelli che
281 // necessitano dell'uso della ricorsione
282 if (j3min==j3max) // big_if
283 v_w3jm[j3max-j3min]=wdown0(j1,j2,m1,m2); // Unico termine da calcolare
284 else {
285 // Si inizializza la ricorsione
286 v_w3jm[j3max-j3min]=wdown0(j1,j2,m1,m2);
287 v_w3jm[j3max-1-j3min]=-(dr(j1,j2,j1+j2,m1,m2,-m1-m2)/( (j1+j2+1)*cr(j1,j2,j1+j2,m1,m2,-m1-m2) ))*v_w3jm[j3max-j3min];
288
289
290 // Ciclo per il calcolo ricorsivo
291 while(j3-2>=j3min){ // down_do
292
293 //Primo coeff della ricorsione
294 double cd1=dr(j1,j2,j3-1,m1,m2,-m1-m2)/(j3*cr(j1,j2,j3-1,m1,m2,-m1-m2));
295 double cd2=((j3-1)*cr(j1,j2,j3,m1,m2,-m1-m2))/(j3*cr(j1,j2,j3-1,m1,m2,-m1-m2));
296 //Ricorsione
297 v_w3jm[j3-2-j3min]=-cd1*v_w3jm[j3-1-j3min]-cd2*v_w3jm[j3-j3min];
298
299 //Aggiorno gli indici
300 --j3;
301 } //END DO down_do
302
303 // Inizializzo gli indici per la upward recursion
304 j3=j3min;
305
306 // Calcolo del primo termine di wigner dal basso
307 v_w3jm[j3-j3min]=wup0(j1,j2,m1,m2);
308
309 // Calcolo del secondo termine di wigner dal basso
310 // Pongo anche la condizione sul coefficienti nel caso ci sia signolarita'
311 double cu3 = (j3min==0) ? 0 : (dr(j1,j2,j3,m1,m2,-m1-m2)/(j3*cr(j1,j2,j3+1,m1,m2,-m1-m2)));
312
313 double w3jm_temp=-cu3*v_w3jm[j3-j3min];
314
315 // If legato alla monotonia della successione
316 if (fabs(w3jm_temp)>fabs(v_w3jm[j3min-j3min])) { //up_if
317
318 // Aggiorno gli indici e metto nell'array il secondo valore
319 // in questo modo sono pronto per iniziale la upward recursion
320 // a tre termini
321 ++j3;
322 v_w3jm[j3-j3min]=w3jm_temp;
323 while(1) { // up_do: DO
324 //Aggiorno gli indici
325 ++j3;
326
327 if (j3-1==j3max) break;
328
329 // IF ((INT(-m1)==-1).AND.(INT(j1)==1).AND.(INT(m2)==1).AND.(INT(j2)==2)) THEN
330 // WRITE(*,*) "j3-1,cr1,cr2",j3-1,cr(j1,j2,j3,m1,m2,-m1-m2),cr(j1,j2,j3,m1,m2,-m1-m2)
331 // END IF
332
333 //Primo e secondo coeff della ricorsione
334 double cu1=dr(j1,j2,j3-1,m1,m2,-m1-m2)/((j3-1)*cr(j1,j2,j3,m1,m2,-m1-m2));
335 double cu2=(j3*cr(j1,j2,j3-1,m1,m2,-m1-m2))/((j3-1)*cr(j1,j2,j3,m1,m2,-m1-m2));
336
337 //Assegnazione temporanea della ricorsione
338 double w3jm_temp=-cu1*v_w3jm[j3-1-j3min]-cu2*v_w3jm[j3-2-j3min];
339
340 if ((fabs(w3jm_temp)<fabs(v_w3jm[j3-1-j3min])) || ((j3-1)==j3max) ) break; // Cond. di uscita
341
342 v_w3jm[j3-j3min]=w3jm_temp; //Assegno perche' e' ok
343 } // END DO up_do
344 } //END IF up_if
345 } // big_if
346 } // END SUBROUTINE wigner3jm
347
348 // calcolo una serie di coefficienti di Gaunt per una data quadrupletta di indici; vec_trans.f90
349 void gaunt_cz(int m, int n, int mu, int nu, int qmaxa, double *v_aq, int *error) {
350 // TODO zrušit error a dát místo něj
351 if (error) *error = 0;
352 if (abs(m) > n || abs(mu) > nu) { // error_if
353 if (error) *error=1;
354 assert(0);
355 return;
356 }
357
358 // calcolo i bounds dei vettori di wigner
359 int pmin = j3minimum(n,nu,m,mu);
360 int pmax = n+nu;
361 int pmin0 = j3minimum(n,nu,0,0);
362 // Alloco i vettori di wigner e li calcolo
363 double *v_w3jm = calloc(pmax-pmin+1, sizeof(double));
364 double *v_w3jm0 = calloc(pmax-pmin+1, sizeof(double));
365 assert(v_w3jm);
366 assert(v_w3jm0);
367 //ALLOCATE(v_w3jm(pmin:pmax),v_w3jm0(pmin0:pmax),STAT=stat_a)
368 wigner3jm(n,nu,m,mu,pmin,pmax,v_w3jm); // FIXME error handling, non-void return values
369 wigner3jm(n,nu,0,0,pmin0,pmax,v_w3jm0);
370
371 // Entro nel ciclo per il calcolo dei coefficienti di gaunt
372 for (int q = 0; q <= qmaxa; ++q) { //gaunt_do: DO q=0,qmaxa
373
374 // Calcolo dell'indice p, sia reale che intero
375 int p=(n+nu)-2*q;
376 //pr=REAL(p,dbl) // use only integer p
377
378 //Calcolo del fattoriale
379 double logw = .5 * (lnf(n+m)+lnf(nu+mu)+lnf(p-m-mu) - \
380 lnf(n-m)-lnf(nu-mu)-lnf(p+m+mu));
381 double fac = exp(logw);
382
383 // Calcolo del coefficiente di gaunt
384 v_aq[q]=min1pow(m+mu)*(2*p+1)*fac*v_w3jm[p]*v_w3jm0[p];
385
386 } // END DO gaunt_do
387
388 // Disalloco i vettori di wigner a lavoro finito
389 free(v_w3jm);
390 free(v_w3jm0);
391 } // END SUBROUTINE gaunt_cz
392
393
394 // gaunt_xu from vec_trans.f90
395 // FIXME set some sensible return value
396 void gaunt_xu(int m, int n, int mu, int nu, int qmax, double *v_aq, int *error) {
397
398 double alphap1;
399 *error = 0;
400 int v_zero[qmax];
401 for (int i = 0; i < qmax; i++) v_zero[i] = 1;
402 double v_aq_cz[qmax];
403 for (int i = 0; i < qmax; i++) v_aq_cz[i] = 0.;
404 int qi = 0;
405
406 if(abs(m)>n || abs(mu)>nu) {
407 *error = 1;
408 fprintf(stderr, "invalid values for m, n, mu or nu\n");
409 return; // FIXME vyřešit chyby
410 }
411
412 switch(qmax) { //qmax_case
413 case 0:
414 v_aq[0] = f_a0(m,n,mu,nu);
415 break;
416 case 1:
417 v_aq[0] = f_a0(m,n,mu,nu);
418 v_aq[1] = v_aq[0] * f_a1norm(m,n,mu,nu);
419 // controllo gli zeri
420 if (fabs(v_aq[1]/v_aq[0]) < ZERO_THRESHOLD) {
421 v_aq[1] = 0.;
422 v_zero[1] = 0;
423 }
424 break;
425 case 2:
426 v_aq[0] = f_a0(m,n,mu,nu);
427
428 v_aq[1] = v_aq[0] * f_a1norm(m,n,mu,nu);
429 // controllo gli zeri
430 if (fabs(v_aq[1]/v_aq[0]) < ZERO_THRESHOLD) {
431 v_aq[1] = 0.;
432 v_zero[1] = 0;
433 }
434
435 v_aq[2] = v_aq[0] * f_a2normr(m,n,mu,nu,v_aq[1]/v_aq[0]);
436 // controllo gli zeri
437 if (fabs(v_aq[2]/v_aq[0]) < ZERO_THRESHOLD) {
438 v_aq[2] = 0.;
439 v_zero[2] = 0;
440 }
441 break;
442 default:
443 if (m == 0 && mu == 0) { // big_if
444 v_aq[0] = f_a0(m,n,mu,nu);
445
446 // backward recursion
447 for (int q = 1; q <= qmax; ++q) { // uno_b_do
448 int p = n + nu - 2*q;
449 double c0 = f_alpha(n,nu,p+1);
450 double c1 = f_alpha(n,nu,p+2);
451 v_aq[q] = (c1/c0) * v_aq[q-1];
452
453 // Vedo se il q-esimo valore e' zero
454 if (v_zero[q-1] == 1) {// v_zero_if_1
455 if(fabs(v_aq[q]/v_aq[q-1]) < ZERO_THRESHOLD) { // zg_if_1
456 v_aq[q] = 0.;
457 v_zero[q] = 0;
458 }
459 } else if (v_zero[q-1]==0 && v_zero[q-2]) {
460 if(fabs(v_aq[q]/v_aq[q-2]) < ZERO_THRESHOLD) {// zg_if1_1
461 v_aq[q] = 0.;
462 v_zero[q] = 0;
463 }
464 } //v_zero_if_1
465 } //uno_b_do
466 } else if (mu == m && nu == n) {
467 v_aq[0] = f_a0(m,n,mu,nu);
468
469 // backward recursion
470 for (int q = 1; q <= qmax; ++q) { // due_b_do
471 // calculate pre-coefficients
472 int p = n + nu - 2*q;
473 int p1 = p - m - mu;
474 int p2 = p + m + mu;
475
476 // calculate recursion coefficients
477 double c0 = (p+2) * (p1+1) * f_alpha(n,nu,p+1);
478 double c1 = (p+1) * (p2+2) * f_alpha(n,nu,p+2);
479
480 //recursion
481 v_aq[q] = (c1/c0) * v_aq[q-1];
482
483 // Vedo se il q-esimo valore e' zero
484 if (v_zero[q-1] == 1) {// v_zero_if_2
485 if(fabs(v_aq[q]/v_aq[q-1]) < ZERO_THRESHOLD) { // zg_if_2
486 v_aq[q] = 0.;
487 v_zero[q] = 0;
488 }
489 } else if (v_zero[q-1]==0 && v_zero[q-2]) {
490 if(fabs(v_aq[q]/v_aq[q-2]) < ZERO_THRESHOLD) {// zg_if1_2
491 v_aq[q] = 0.;
492 v_zero[q] = 0;
493 }
494 } //v_zero_if_2
495 } // due_b_do
496 } else if (mu == -m) {
497 // Primo valore per la backward recursion
498 v_aq[0] = f_a0(m,n,mu,nu);
499 v_aq[1] = f_a1norm(m,n,mu,nu)*v_aq[0];
500
501 // Controllo gli zeri
502 if (fabs(v_aq[1]/v_aq[0]) < ZERO_THRESHOLD) { //zg_if_3_0
503 v_aq[1] = 0.;
504 v_zero[1] = 0;
505 } //zg_if_3_0
506
507 // backward recursion
508 for (int q = 2; q <= qmax; ++q) { // tre_b_do
509 // calculate pre-coefficient
510 int p = n + nu - 2*q;
511
512 // calculate recursion coefficients
513 double c0 = f_alpha(n, nu, p+1);
514 double c1 = 4*isq(m) + f_alpha(n,nu,p+2) + f_alpha(n,nu,p+3);
515 double c2 = - f_alpha(n, nu, p+4);
516
517 // recursion
518 v_aq[q] = (c1/c0)*v_aq[q-1] + (c2/c0)*v_aq[q-2];
519
520 // Vedo se il q-esimo valore e' zero
521 if (v_zero[q-1] == 1) {// v_zero_if_3
522 if(fabs(v_aq[q]/v_aq[q-1]) < ZERO_THRESHOLD) { // zg_if_3
523 v_aq[q] = 0.;
524 v_zero[q] = 0;
525 }
526 } else if (v_zero[q-1]==0 && v_zero[q-2]) {
527 if(fabs(v_aq[q]/v_aq[q-2]) < ZERO_THRESHOLD) {// zg_if1_3
528 v_aq[q] = 0.;
529 v_zero[q] = 0;
530 }
531 } //v_zero_if_3
532 } // tre_b_do
533
534 // forward recursion
535 // Primo valore per la forward recursion,errore relativo e suo swap
536 double aq_fwd=f_aqmax(m,n,mu,nu,qmax);
537 double res=fabs(aq_fwd-v_aq[qmax])/fabs(aq_fwd);
538
539 //Se non ho precisione, sostituisco i valori
540 if (res>BF_PREC) { //tre_f_if
541 v_aq[qmax]=aq_fwd;
542 int qi=1;
543 int zeroswitch = 0; // black magic (gevero's "switch")
544 //Entro nel ciclo della sostituzione valori
545 for( int q=qmax-1;q>=0;--q) { // tre_f_do
546 switch(qmax-q) {// tre_q_case // FIXME switch -> if/else
547 case 1: {// q==qmax-1
548 //Calcolo v_aq[qmax-1]
549 int p=n+nu-2*(q+2);
550 double c1=4*isq(m)+f_alpha(n,nu,p+2)+f_alpha(n,nu,p+3);
551 double c2=-f_alpha(n,nu,p+4);
552 double aq_fwd=-(c1/c2)*v_aq[qmax];
553
554 switch(v_zero[q]) { // z_3_1_case
555 case 0:
556 v_aq[q] = 0.;
557 break;
558 case 1:
559 res=fabs(aq_fwd-v_aq[q])/fabs(aq_fwd);
560 break;
561 default:
562 assert(0);
563 }
564 }
565 break;
566 default: { //Per tutti gli altri q
567 //Calcolo v_aq[qmax-1]
568 int p=n+nu-2*(q+2);
569 double c0=f_alpha(n,nu,p+1);
570 double c1=4*isq(m)+f_alpha(n,nu,p+2)+f_alpha(n,nu,p+3);
571 double c2=-f_alpha(n,nu,p+4);
572 aq_fwd=-(c1/c2)*v_aq[q+1]+(c0/c2)*v_aq[q+2];
573
574 switch(v_zero[q]){ // z_3_2_case//FIXME switch -> if/else
575 case 0:
576 v_aq[q] = 0.;
577 break;
578 case 1: //Il valore precedente e' zero
579 res=fabs(aq_fwd-v_aq[q])/fabs(aq_fwd);
580 break;
581 default:
582 assert(0);
583 }
584 }
585 } //tre_q_case
586 //Adesso se la precisione e' raggiunta esco dal ciclo,
587 //se no sostituisco e rimango
588 if (res<BF_PREC || q==0 || fabs(aq_fwd) < fabs(v_aq[q+1]))
589 break; //tre_f_do
590 //Sono nel ciclo, allora sostituisco eaggiorno indice e residuo
591 v_aq[q]=aq_fwd;
592 qi=q;
593 assert(q); // assert níže přesunut sem
594 } // tre_f_do
595 // Check sul ciclo di sostituzione
596 //assert(q);
597 /*
598 error_if1: IF (q==0) THEN
599 WRITE(*,*)
600 WRITE(*,*) "Si e' verificato un errore nella subroutine gaunt_xu:"
601 WRITE(*,*) "la precisione richiesta per i coefficienti di Gaunt nella backward"
602 WRITE(*,*) "e forward recursion non e' stata raggiunta"
603 WRITE(*,*)
604 error=1
605 RETURN
606 END IF error_if1
607 */
608 } // tre_f_if
609 } else { // caso generale (4)
610 // backward
611 // Calcolo direttamente i primi tre valori della ricorsione
612 v_aq[0]=f_a0(m,n,mu,nu);
613 v_aq[1]=v_aq[0]*f_a1norm(m,n,mu,nu);
614
615 // vedo se il secondo valore e' zero
616 if (fabs(v_aq[1]/v_aq[0]) < ZERO_THRESHOLD) { // zg1_if
617 v_aq[1] = 0.;
618 v_zero[1] = 0;
619 }
620
621 //...........................................................
622 //Calcolo il terzo valore della ricorsione in funzione di Ap4
623 //...........................................................
624 //Inizializzo i valori comuni per i coefficienti
625 int p=n+nu-2*(2);
626 int p1=p-m-mu;
627 int p2=p+m+mu;
628 alphap1=f_alpha(n,nu,p+1);
629 double alphap2=f_alpha(n,nu,p+2);
630 double Ap2=f_Ap(m,n,mu,nu,p+2);
631 double Ap3=f_Ap(m,n,mu,nu,p+3);
632 double Ap4=f_Ap(m,n,mu,nu,p+4);
633
634 //Con questo if decido se mi serve la ricorsione a 3 o 4 termini
635 if (Ap4==0) { //Ap4_2_if
636 //Calcolo i restanti valori preliminari
637 double Ap5=f_Ap(m,n,mu,nu,p+5);
638 double Ap6=f_Ap(m,n,mu,nu,p+6);
639 double alphap5=f_alpha(n,nu,p+5);
640 double alphap6=f_alpha(n,nu,p+6);
641
642 //Calcolo i coefficienti per la ricorsione ma non c3 perche' qui e solo qui non mi serve
643 double c0=(p+2.)*(p+3.)*(p+5.)*(p1+1.)*(p1+2.)*(p1+4.)*Ap6*alphap1;
644 double c1=(p+5.)*(p1+4.)*Ap6*(Ap2*Ap3 + (p+1.)*(p+3.)*(p1+2.)*(p2+2.)*alphap2);
645 double c2=(p+2.)*(p2+3.)*Ap2*(Ap5*Ap6 + (p+4.)*(p+6.)*(p1+5.)*(p2+5.)*alphap5);
646
647 //Calcolo il mio coefficiente
648 v_aq[2]=(c1/c0)*v_aq[1]+(c2/c0)*v_aq[0];
649
650 //Assegno l'indice segnaposto per Ap4=0
651 // q4=2 FIXME UNUSED
652 } else {
653 //Calcolo i restanti valori preliminari
654 double alphap3=f_alpha(n,nu,p+3);
655 double alphap4=f_alpha(n,nu,p+4);
656
657 //Calcolo coefficienti ricorsione
658 double c0=(p+2.)*(p+3.)*(p1+1.)*(p1+2.)*Ap4*alphap1;
659 double c1=Ap2*Ap3*Ap4+(p+1.)*(p+3.)*(p1+2.)*(p2+2.)*Ap4*alphap2+ \
660 (p+2.)*(p+4.)*(p1+3.)*(p2+3.)*Ap2*alphap3;
661 double c2=-(p+2.)*(p+3.)*(p2+3.)*(p2+4.)*Ap2*alphap4;
662
663 //Calcolo il mio coefficiente
664 v_aq[2]=(c1/c0)*v_aq[1]+(c2/c0)*v_aq[0];
665 } // Ap4_2_if
666
667 //Vedo se il terzo valore e' zero
668 if (v_zero[1]==1) { // v_zero_if1
669 if (fabs(v_aq[2]/v_aq[1])< ZERO_THRESHOLD) { //zg2_if
670 v_aq[2]=0;
671 v_zero[2]=0;
672 }
673 } else if (v_zero[1]==0) {
674 if (fabs(v_aq[2]/v_aq[0])<ZERO_THRESHOLD) { //zg2_if1:
675 v_aq[2]=0;
676 v_zero[2]=0;
677 }
678 } // v_zero_if1
679
680
681 //...........................................................
682 //Calcolo i restanti valori nel loop
683 //...........................................................
684 for (int q = 3; q <= qmax; q++ ) { //gen_bwd_do: DO q=3,qmax
685
686 //Inizializzo i valori comuni per i coefficienti
687 int p=n+nu-2*(q);
688 int p1=p-m-mu;
689 int p2=p+m+mu;
690 alphap1=f_alpha(n,nu,p+1);
691 double alphap2=f_alpha(n,nu,p+2);
692 double Ap2=f_Ap(m,n,mu,nu,p+2);
693 double Ap3=f_Ap(m,n,mu,nu,p+3);
694 double Ap4=f_Ap(m,n,mu,nu,p+4);
695
696 //Con questo if decido se mi serve la ricorsione a 3 o 4 termini
697 if (Ap4==0) { // Ap4_bwd_if: IF (Ap4==0) THEN
698
699 //Calcolo i restanti valori preliminari
700 double Ap5=f_Ap(m,n,mu,nu,p+5);
701 double Ap6=f_Ap(m,n,mu,nu,p+6);
702 double alphap5=f_alpha(n,nu,p+5);
703 double alphap6=f_alpha(n,nu,p+6);
704
705 //Calcolo i coefficienti per la ricorsione ma non c3 perche' qui e solo qui non mi serve
706 double c0=(p+2.)*(p+3.)*(p+5.)*(p1+1.)*(p1+2.)*(p1+4.)*Ap6*alphap1;
707 double c1=(p+5.)*(p1+4.)*Ap6*(Ap2*Ap3 + (p+1.)*(p+3.)*(p1+2.)*(p2+2.)*alphap2);
708 double c2=(p+2.)*(p2+3.)*Ap2*(Ap5*Ap6 + (p+4.)*(p+6.)*(p1+5.)*(p2+5.)*alphap5);
709 double c3=-(p+2.)*(p+4.)*(p+5.)*(p2+3.)*(p2+5.)*(p2+6.)*Ap2*alphap6;
710
711 //Calcolo il mio coefficiente
712 v_aq[q]=(c1/c0)*v_aq[q-1]+(c2/c0)*v_aq[q-2]+(c3/c0)*v_aq[q-3];
713
714 //Assegno l'indice segnaposto per Ap4=0
715 //q4=q // FIXME nepoužitá proměnná
716 } else {
717 //Calcolo i restanti valori preliminari
718 double alphap3=f_alpha(n,nu,p+3);
719 double alphap4=f_alpha(n,nu,p+4);
720
721 //Calcolo coefficienti ricorsione
722 double c0=(p+2.)*(p+3.)*(p1+1.)*(p1+2.)*Ap4*alphap1;
723 double c1=Ap2*Ap3*Ap4+(p+1.)*(p+3.)*(p1+2.)*(p2+2.)*Ap4*alphap2+ \
724 (p+2.)*(p+4.)*(p1+3.)*(p2+3.)*Ap2*alphap3;
725 double c2=-(p+2.)*(p+3.)*(p2+3.)*(p2+4.)*Ap2*alphap4;
726
727 //Calcolo il mio coefficiente
728 v_aq[q]=(c1/c0)*v_aq[q-1]+(c2/c0)*v_aq[q-2];
729 } // END IF Ap4_bwd_if
730
731 //Vedo se il q-esimo valore e' zero
732 if(v_zero[q-1]==1) { //v_zero_ifq: IF (v_zero[q-1]==1) THEN
733 if(fabs(v_aq[q]/v_aq[q-1])<ZERO_THRESHOLD) {//zgq_if
734 v_aq[q]=0.;
735 v_zero[q]=0;
736 }
737 } else if ((v_zero[q-1]==0) && (v_zero[q-2] !=0)) {
738 if (fabs(v_aq[q]/v_aq[q-2])<ZERO_THRESHOLD) {//zgq_if1:
739 v_aq[q]=0.;
740 v_zero[q]=0;
741 } // zgq_if1
742 } // v_zero_ifq
743
744 } //gen_bwd_do
745
746 //---------------------------------------------------------------------------------
747 //FORWARD
748 //---------------------------------------------------------------------------------
749
750 //Calcolo pmin,Apmin e la mia variabile logica
751 int pmin=n+nu-2*(qmax);
752 int Apmin=f_Ap(m,n,mu,nu,pmin);
753 int test=(((Apmin)==0) &&
754 (((pmin)==(m+mu+1)) || ((pmin)==(-m-mu+1))));
755
756 //........................................................
757 //Se la mia variabile logica e' vera, Faccio il mio conto
758 //........................................................
759 if(test) { //Apmin_if: if (test) THEN
760
761 //Il valore per qmax allora e' zero
762 v_aq[qmax]=0;
763
764 //Calcolo il secondo valore, e se la precisione e' raggiunta esco
765 double aq_fwd=f_aqmax_1(m,n,mu,nu,qmax);
766 double res=fabs(aq_fwd-v_aq[qmax-1])/fabs(aq_fwd);
767 if (res<BF_PREC)
768 return;
769
770 //Assegno il secondo valore e faccio il ciclo
771 v_aq[qmax-1]=aq_fwd;
772 qi=1; //FIXME nepoužitá proměnná
773
774 for (int q = qmax; q >= 2; --q) { //Apmin_do: DO q=qmax,2,-1
775
776 //Calcolo pre-coefficienti
777 int p=n+nu-2*(q);
778 int p1=p-m-mu;
779 int p2=p+m+mu;
780 alphap1=f_alpha(n,nu,p+1);
781 double alphap2=f_alpha(n,nu,p+2);
782 double alphap3=f_alpha(n,nu,p+3);
783 double alphap4=f_alpha(n,nu,p+4);
784 double Ap2=f_Ap(m,n,mu,nu,p+2);
785 double Ap3=f_Ap(m,n,mu,nu,p+3);
786 double Ap4=f_Ap(m,n,mu,nu,p+4);
787
788 //Calcolo coefficienti ricorsione
789 double c0=(p+2.)*(p+3.)*(p1+1.)*(p1+2.)*Ap4*alphap1;
790 double c1=Ap2*Ap3*Ap4+(p+1.)*(p+3.)*(p1+2.)*(p2+2.)*Ap4*alphap2+
791 (p+2.)*(p+4.)*(p1+3.)*(p2+3.)*Ap2*alphap3;
792 double c2=-(p+2.)*(p+3.)*(p2+3.)*(p2+4.)*Ap2*alphap4;
793
794 //Ricorsione e residuo
795 aq_fwd=-(c1/c2)*v_aq[q-1]+(c0/c2)*v_aq[q];
796 res=fabs(aq_fwd-v_aq[q-2])/fabs(aq_fwd);
797
798 if (res<BF_PREC) return;
799
800 v_aq[q-2]=aq_fwd;
801 qi=q-2;
802
803 } // END DO Apmin_do
804
805 // Check sul ciclo di sostituzione
806 assert (qi); // Apmin_error_if1
807 /*{
808 WRITE(*,*)
809 WRITE(*,*) "Si e' verificato un errore nella subroutine gaunt_xu, caso generale, Apmin=0:"
810 WRITE(*,*) "la precisione richiesta per i coefficienti di Gaunt nella backward"
811 WRITE(*,*) "e forward recursion non e' stata raggiunta"
812 WRITE(*,*)
813 error=1
814 RETURN
815 }*/ // Apmin_error_if1
816
817 //Esco dalla subroutine gaunt_xu
818 return;
819
820 } // Apmin_if
821
822 //..........................................................................
823 //CASO GENERALE PER LA FORWARD RECURRENCE
824 //..........................................................................
825
826 //Primo valore per la forward recursion,errore relativo e suo swap
827 double aq_fwd=f_aqmax(m,n,mu,nu,qmax);
828 double res=fabs(aq_fwd-v_aq[qmax])/fabs(aq_fwd);
829 qi=1;
830
831 if (res>BF_PREC) { //gen_f_if
832 //Se non ho precisione, sostituisco i valori
833
834 v_aq[qmax]=aq_fwd;
835
836 qi=qmax-1;
837
838 //Entro nel ciclo della sostituzione valori
839 while(1) { // gen_f_do: DO
840
841 switch(qmax-qi) {//gen_q_case:SELECT CASE (qmax-qi)
842
843 //$$$$$$$$$$$$$$$$
844 case 1: { //q=qmax-1
845 //$$$$$$$$$$$$$$$$
846 //Calcolo Ap4 per qi+2 per vedere quale schema usare
847 int p=n+nu-2*(qi+2);
848 double Ap4=f_Ap(m,n,mu,nu,p+4);
849
850 //Scelgo la ricorrenza a seconda del valore di Ap4
851 if (Ap4==0) { // Ap4_q1_if
852
853 //Calcolo aq secondo la ricorrenza a 4 termini: uso qi+3 perche' il termine piu' alto e'
854 //maggiore di 3 unita' rispetto a qi, pur essendo nullo e non comparendo nella ricorsione
855 int p=n+nu-2*(qi+3);
856 int p1=p-m-mu;
857 int p2=p+m+mu;
858 double alphap5=f_alpha(n,nu,p+5);
859 double alphap6=f_alpha(n,nu,p+6);
860 double Ap2=f_Ap(m,n,mu,nu,p+2);
861 double Ap5=f_Ap(m,n,mu,nu,p+5);
862 double Ap6=f_Ap(m,n,mu,nu,p+6);
863 double c2=(p+2.)*(p2+3.)*Ap2*(Ap5*Ap6 + (p+4.)*(p+6.)*(p1+5.)*(p2+5.)*alphap5);
864 double c3=-(p+2.)*(p+4.)*(p+5.)*(p2+3.)*(p2+5.)*(p2+6.)*Ap2*alphap6;
865 aq_fwd=-(c2/c3)*v_aq[qi+1];
866
867 //A seconda che il mio valore sia 0 o meno confronto i valori
868 switch(v_zero[qi]){ //zAp41_case:SELECT CASE (v_zero[qi])
869 case 0:
870 v_aq[qi]=0;
871 break;
872 case 1:
873 res=fabs(aq_fwd-v_aq[qi])/fabs(aq_fwd);
874 if (res<BF_PREC) {
875 //EXIT gen_f_do
876 //měli bychom breaknout smyčku, ale za
877 //ní už nic „smysluplného“ není
878 assert(qi);
879 return;
880 }
881 break;
882 default:
883 assert(0);
884 } // END SELECT zAp41_case
885
886 //Qui calcolo il valore successivo dopo aver aggiornato qi:
887 //Se v_aq[qi]=0 allora non chiamo cruzan, se no lo chamo e
888 //tengo un solo valore
889 qi=qi-1;
890
891 switch(v_zero[qi]) {//zcz1_case:SELECT CASE (v_zero[qi])
892 case 0:
893 v_aq[qi]=0;
894 qi=qi-1;
895 continue; //CYCLE gen_f_do
896 break;
897 case 1:
898 gaunt_cz(m,n,mu,nu,qmax,&(v_aq_cz[qi]),error); // FIXME implementace gaunt_cz
899 aq_fwd=(v_aq_cz[qi]);
900 res=fabs(aq_fwd-v_aq[qi])/fabs(aq_fwd);
901 break;
902 default:
903 assert(0);
904 }
905
906 //-----------------
907 } else { //Qui Ap4/=0
908 //-----------------
909 //Calcolo aq
910 int p=n+nu-2*(qi+2);
911 int p1=p-m-mu;
912 int p2=p+m+mu;
913 double alphap2=f_alpha(n,nu,p+2);
914 double alphap3=f_alpha(n,nu,p+3);
915 double alphap4=f_alpha(n,nu,p+4);
916 double Ap2=f_Ap(m,n,mu,nu,p+2);
917 double Ap3=f_Ap(m,n,mu,nu,p+3);
918 double Ap4=f_Ap(m,n,mu,nu,p+4);
919 double c1=Ap2*Ap3*Ap4+(p+1.)*(p+3.)*(p1+2.)*(p2+2.)*Ap4*alphap2+
920 (p+2.)*(p+4.)*(p1+3.)*(p2+3.)*Ap2*alphap3;
921 double c2=-(p+2.)*(p+3.)*(p2+3.)*(p2+4.)*Ap2*alphap4;
922 aq_fwd=-(c1/c2)*v_aq[qi+1]; //E' qui che lo calcolo
923
924 //A seconda che il mio valore sia 0 o meno confronto i valori
925 switch(v_zero[qi]) { // zAp4d1_case:SELECT CASE (v_zero[qi])
926 case 0:
927 v_aq[qi]=0;
928 qi=qi-1;
929 continue; // gen_f_do
930 break;
931 case 1:
932 res=fabs(aq_fwd-v_aq[qi])/fabs(aq_fwd);
933 break;
934 default:
935 assert(0);
936 }
937
938 } // Ap4_q1_if
939 } break;
940
941 //$$$$$$$$$$$$$$$$
942 case 2: {//CASE(2) gen_q_case //q=qmax-2
943 //$$$$$$$$$$$$$$$$
944
945
946
947
948
949
950 //Calcolo Ap4 per qi+2 per vedere quale schema usare
951 int p=n+nu-2*(qi+2);
952 double Ap4=f_Ap(m,n,mu,nu,p+4);
953
954 //Scelgo la ricorrenza a seconda del valore di Ap4
955 if (Ap4==0) { // Ap4_q2_if
956 //Calcolo aq secondo la ricorrenza a 4 termini: uso qi+3 perche' il termine piu' alto e'
957 //maggiore di 3 unita' rispetto a qi, pur essendo nullo e non comparendo nella ricorsione
958 int p=n+nu-2*(qi+3);
959 int p1=p-m-mu;
960 int p2=p+m+mu;
961 double alphap2=f_alpha(n,nu,p+2);
962 double alphap5=f_alpha(n,nu,p+5);
963 double alphap6=f_alpha(n,nu,p+6);
964 double Ap2=f_Ap(m,n,mu,nu,p+2);
965 double Ap3=f_Ap(m,n,mu,nu,p+3);
966 double Ap5=f_Ap(m,n,mu,nu,p+5);
967 double Ap6=f_Ap(m,n,mu,nu,p+6);
968 double c1=(p+5.)*(p1+4.)*Ap6*(Ap2*Ap3 + (p+1.)*(p+3.)*(p1+2.)*(p2+2.)*alphap2);
969 double c2=(p+2.)*(p2+3.)*Ap2*(Ap5*Ap6 + (p+4.)*(p+6.)*(p1+5.)*(p2+5.)*alphap5);
970 double c3=-(p+2.)*(p+4.)*(p+5.)*(p2+3.)*(p2+5.)*(p2+6.)*Ap2*alphap6;
971 aq_fwd=-(c1/c3)*v_aq[qi+2] -(c2/c3)*v_aq[qi+1];
972
973 //A seconda che il mio valore sia 0 o meno confronto i valori
974 switch(v_zero[qi]) { // zAp42_case
975 case 0:
976 v_aq[qi]=0;
977 break;
978 case 1:
979 res=fabs(aq_fwd-v_aq[qi])/fabs(aq_fwd);
980 if (res<BF_PREC) { // EXIT gen_f_do
981 assert(qi);
982 return;
983 }
984 break;
985 default:
986 assert(0);
987 } // zAp42_case
988
989 //Qui calcolo il valore successivo dopo aver aggiornato qi:
990 //Se v_aq[qi]=0 allora non chiamo cruzan, se no lo chamo e
991 //tengo un solo valore
992 qi=qi-1;
993
994 switch (v_zero[qi]) {//zcz2_case:SELECT CASE (v_zero[qi])
995 case 0:
996 v_aq[qi]=0;
997 qi=qi-1;
998 continue; // gen_f_do
999 break;
1000 case 1:
1001 //FIXME gaunt_cz !!!!!!!!!!!!!!!!!!
1002 gaunt_cz(m,n,mu,nu,qmax,&(v_aq_cz[qi]),error);
1003 aq_fwd=v_aq_cz[qi];
1004 res=fabs(aq_fwd-v_aq[qi])/fabs(aq_fwd);
1005 break;
1006 default:
1007 assert(0);
1008 }
1009 } else { //Qui Ap4!=0
1010 //Calcolo aq
1011 int p=n+nu-2*(qi+2);
1012 int p1=p-m-mu;
1013 int p2=p+m+mu;
1014 double alphap2=f_alpha(n,nu,p+2);
1015 double alphap3=f_alpha(n,nu,p+3);
1016 double alphap4=f_alpha(n,nu,p+4);
1017 double Ap2=f_Ap(m,n,mu,nu,p+2);
1018 double Ap3=f_Ap(m,n,mu,nu,p+3);
1019 double Ap4=f_Ap(m,n,mu,nu,p+4);
1020 double c0=(p+2.)*(p+3.)*(p1+1.)*(p1+2.)*Ap4*alphap1;
1021 double c1=Ap2*Ap3*Ap4+(p+1.)*(p+3.)*(p1+2.)*(p2+2.)*Ap4*alphap2+
1022 (p+2.)*(p+4.)*(p1+3.)*(p2+3.)*Ap2*alphap3;
1023 double c2=-(p+2.)*(p+3.)*(p2+3.)*(p2+4.)*Ap2*alphap4;
1024 aq_fwd=(c0/c2)*v_aq[qi+2]-(c1/c2)*v_aq[qi+1]; //E' qui che lo calcolo
1025
1026 //A seconda che il mio valore sia 0 o meno confronto i valori
1027 switch(v_zero[qi]) { //zAp4d2_case:SELECT CASE (v_zero[qi])
1028 case 0:
1029 v_aq[qi]=0;
1030 qi=qi-1;
1031 continue; // gen_f_do
1032 break;
1033 case 1:
1034 res=fabs(aq_fwd-v_aq[qi])/fabs(aq_fwd);
1035 break;
1036 default:
1037 assert(0);
1038 }
1039
1040 } // Ap4_q2_if
1041 } break;
1042
1043
1044 //$$$$$$$$$$$$$$$$$$$$$$
1045 default: { //CASE DEFAULT gen_q_case //Per tutti gli altri q
1046 //$$$$$$$$$$$$$$$$$$$$$$
1047
1048 //Calcolo Ap4 per qi+2 per vedere quale schema usare
1049 int p=n+nu-2*(qi+2);
1050 double Ap4=f_Ap(m,n,mu,nu,p+4);
1051
1052 //Scelgo la ricorrenza a seconda del valore di Ap4
1053 if (Ap4==0) { // Ap4_qq_if
1054
1055 //Calcolo aq secondo la ricorrenza a 4 termini: uso qi+3 perche' il termine piu' alto e'
1056 //maggiore di 3 unita' rispetto a qi, pur essendo nullo e non comparendo nella ricorsione
1057 int p=n+nu-2*(qi+3);
1058 int p1=p-m-mu;
1059 int p2=p+m+mu;
1060 double alphap2=f_alpha(n,nu,p+2);
1061 double alphap5=f_alpha(n,nu,p+5);
1062 double alphap6=f_alpha(n,nu,p+6);
1063 double Ap2=f_Ap(m,n,mu,nu,p+2);
1064 double Ap3=f_Ap(m,n,mu,nu,p+3);
1065 double Ap5=f_Ap(m,n,mu,nu,p+5);
1066 double Ap6=f_Ap(m,n,mu,nu,p+6);
1067 double c0=(p+2.)*(p+3.)*(p+5.)*(p1+1.)*(p1+2.)*(p1+4.)*Ap6*alphap1;
1068 double c1=(p+5.)*(p1+4.)*Ap6*(Ap2*Ap3 + (p+1.)*(p+3.)*(p1+2.)*(p2+2.)*alphap2);
1069 double c2=(p+2.)*(p2+3.)*Ap2*(Ap5*Ap6 + (p+4.)*(p+6.)*(p1+5.)*(p2+5.)*alphap5);
1070 double c3=-(p+2.)*(p+4.)*(p+5.)*(p2+3.)*(p2+5.)*(p2+6.)*Ap2*alphap6;
1071 aq_fwd=(c0/c3)*v_aq[qi+3]-(c1/c3)*v_aq[qi+2] -(c2/c3)*v_aq[qi+1];
1072
1073 //A seconda che il mio valore sia 0 o meno confronto i valori
1074 switch(v_zero[qi]) {//zAp4q_case:SELECT CASE (v_zero[qi])
1075 case 0:
1076 v_aq[qi]=0;
1077 break;
1078 case 1:
1079 res=fabs(aq_fwd-v_aq[qi])/fabs(aq_fwd);
1080 if (res<BF_PREC) {// EXIT gen_f_do
1081 assert(qi);
1082 return;
1083 }
1084 break;
1085 default:
1086 assert(0);
1087 }
1088
1089 //Qui calcolo il valore successivo dopo aver aggiornato qi:
1090 //Se v_aq[qi]=0 allora non chiamo cruzan, se no lo chiamo e
1091 //tengo un solo valore.L'if c'e' per non far sballare qi
1092
1093 if (qi>0) { // qi_if
1094
1095 qi=qi-1;
1096
1097 switch(v_zero[qi]) { //zczq_case:SELECT CASE (v_zero[qi])
1098 case 0:
1099 v_aq[qi]=0;
1100 qi=qi-1;
1101 continue; // CYCLE gen_f_do
1102 break;
1103 case 1:
1104 gaunt_cz(m,n,mu,nu,qmax,&(v_aq_cz[qi]),error); // FIXME je potřeba mít v_aq_cz jako pole?
1105 aq_fwd=v_aq_cz[qi];
1106 res=fabs(aq_fwd-v_aq[qi])/fabs(aq_fwd);
1107 break;
1108 default:
1109 assert(0);
1110 }
1111
1112 } // qi_if
1113
1114 //-----------------
1115 } else { //Qui Ap4/=0
1116 //-----------------
1117
1118 //Calcolo aq
1119 int p=n+nu-2*(qi+2);
1120 int p1=p-m-mu;
1121 int p2=p+m+mu;
1122 double alphap2=f_alpha(n,nu,p+2);
1123 double alphap3=f_alpha(n,nu,p+3);
1124 double alphap4=f_alpha(n,nu,p+4);
1125 double Ap2=f_Ap(m,n,mu,nu,p+2);
1126 double Ap3=f_Ap(m,n,mu,nu,p+3);
1127 double Ap4=f_Ap(m,n,mu,nu,p+4);
1128 double c0=(p+2.)*(p+3.)*(p1+1.)*(p1+2.)*Ap4*alphap1;
1129 double c1=Ap2*Ap3*Ap4+(p+1.)*(p+3.)*(p1+2.)*(p2+2.)*Ap4*alphap2+
1130 (p+2.)*(p+4.)*(p1+3.)*(p2+3.)*Ap2*alphap3;
1131 double c2=-(p+2.)*(p+3.)*(p2+3.)*(p2+4.)*Ap2*alphap4;
1132 aq_fwd=(c0/c2)*v_aq[qi+2]-(c1/c2)*v_aq[qi+1]; //E' qui che lo calcolo
1133
1134 //A seconda che il mio valore sia 0 o meno confronto i valori
1135 switch(v_zero[qi]) { //zAp4dq_case:SELECT CASE (v_zero[qi])
1136 case 0:
1137 v_aq[qi]=0;
1138 qi=qi-1;
1139 continue; // gen_f_do
1140 break;
1141 case 1:
1142 res=fabs(aq_fwd-v_aq[qi])/fabs(aq_fwd);
1143 break;
1144 default:
1145 assert(0);
1146 }
1147
1148 } // Ap4_qq_if
1149 } // default
1150
1151 } // END SELECT gen_q_case
1152
1153 //Adesso se la precisione e' raggiunta esco dal ciclo, se no sostituisco e rimango
1154 if ((res<BF_PREC) || (qi==0) || (fabs(aq_fwd)<fabs(v_aq[qi+1]))) // EXIT
1155 break; // gen_f_do
1156
1157 //Sono nel ciclo, allora sostituisco eaggiorno indice e residuo
1158 v_aq[qi]=aq_fwd;
1159 qi=qi-1;
1160
1161 } // END DO gen_f_do
1162
1163 // Check sul ciclo di sostituzione
1164 assert(qi); /* gen_error_if1: if (qi==0) {
1165 WRITE(*,*)
1166 WRITE(*,*) "Si e' verificato un errore nella subroutine gaunt_xu,caso generale:"
1167 WRITE(*,*) "la precisione richiesta per i coefficienti di Gaunt nella backward"
1168 WRITE(*,*) "e forward recursion non e' stata raggiunta"
1169 WRITE(*,*)
1170 error=1
1171 RETURN
1172 } */ // gen_error_if1
1173
1174
1175 } // gen_f_if
1176
1177 } // big_if
1178 } // qmax_case
1179 } // gaunt_xu
1180
1181
1182 #define MIN(x,y) (((x) > (y)) ? (y) : (x))
1183 static inline int q_max(int m, int n, int mu, int nu) {
1184 return MIN(n, MIN(nu,(n+nu-abs(m+mu))/2));
1185 }
1186
1187 /* THIS THING IS TOTALLY WRONG
1188 int gaunt(int m, int n, int mu, int nu, double *v) {
1189 int err = 0;
1190 int qmax = q_max(m,n,mu,nu);
1191 if (!v) v = calloc(qmax+1, sizeof(double));
1192 if (!v) return -1;
1193 gaunt_xu(m, n, mu, nu, qmax, v, &err);
1194 return err;
1195 }
1196 */
1197
1198 #ifdef GAUNTTEST
1199
1200 void __vec_trans_MOD_gaunt_xu(const double *m, const double *n, const double *mu, const double *nu, const int *qmax, double *v_aq, int *err);
1201 int main()
1202 {
1203 int NMAX=20, REPEAT=10;
1204 //int scanned;
1205 int m, n, mu, nu;
1206 for (int r = 0; r < REPEAT; ++r) for (n = 1; n < NMAX; ++n) for (nu = 1 ; nu < NMAX; ++nu)
1207 for (mu = -nu; mu <=nu; ++mu) for (m = -n; m <= n; ++m) {
1208 //while(EOF != (scanned = scanf("%d %d %d %d", &m, &n, &mu, &nu))) {
1209 // if (scanned != 4) continue;
1210 // printf("%d %d %d %d\n", m, n, mu, nu);
1211 double mf = m, nf = n, muf = mu, nuf = nu;
1212 int qmax = q_max(m,n,mu,nu);
1213 double v_aq_c[qmax+1], v_aq_f[qmax+1];
1214 int errc, errf;
1215 __vec_trans_MOD_gaunt_xu(&mf, &nf, &muf, &nuf, &qmax, v_aq_f, &errf);
1216 gaunt_xu(m,n,mu,nu,qmax,v_aq_c, &errc);
1217 // for(int i = 0; i <= qmax; ++i) printf("%f ", v_aq_c[i]);
1218 // puts("(c)");
1219 // for(int i = 0; i <= qmax; ++i) printf("%f ", v_aq_f[i]);
1220 // puts("(f)");
1221 // for(int i = 0; i <= qmax; ++i) printf("%e ", v_aq_f[i] - v_aq_c[i]);
1222 // puts("(f-c)");
1223
1224 }
1225 return 0;
1226 }
1227 #endif //GAUNTTEST
1228
1229 #endif //USE_FORTRAN_GAUNT_XU
QPMS photonic multiple scattering suite