funct_fermi1.f90
Go to the documentation of this file.
1 !
2 ! file contains fermi function routines:
3 !
4 ! function zfermim12 does a rational function fit for the order -1/2 integral
5 ! function zfermi12 does a rational function fit for the order 1/2 integral
6 ! function zfermi1 does a rational function fit for the order 1 integral
7 ! function zfermi32 does a rational function fit for the order 3/2 integral
8 ! function zfermi2 does a rational function fit for the order 2 integral
9 ! function zfermi52 does a rational function fit for the order 5/2 integral
10 ! function zfermi3 does a rational function fit for the order 3 integral
11 !
12 ! function ifermim12 is a rational function fit for the inverse of order -1/2
13 ! function ifermi12 is a rational function fit for the inverse of order 1/2
14 ! function ifermi32 is a rational function fit for the inverse of order 3/2
15 ! function ifermi52 is a rational function fit for the inverse of order 5/2
16 
17 
18 
19 
20 
21  real*8 function zfermim12(x)
22  include 'implno.dek'
23 !
24 ! this routine applies a rational function expansion to get the fermi function
25 ! of order -1/2 evaluated at x. maximum error is 1.23d-12.
26 ! reference: antia apjs 84,101 1993
27 !
28 ! declare
29  integer :: i,m1,k1,m2,k2
30  real*8 :: x,an,a1(12),b1(12),a2(12),b2(12),rn,den,xx
31 
32 ! load the coefficients of the expansion
33  data an,m1,k1,m2,k2 /-0.5d0, 7, 7, 11, 11/
34  data (a1(i),i=1,8)/ 1.71446374704454d7, 3.88148302324068d7, &
35  3.16743385304962d7, 1.14587609192151d7, &
36  1.83696370756153d6, 1.14980998186874d5, &
37  1.98276889924768d3, 1.0d0/
38  data (b1(i),i=1,8)/ 9.67282587452899d6, 2.87386436731785d7, &
39  3.26070130734158d7, 1.77657027846367d7, &
40  4.81648022267831d6, 6.13709569333207d5, &
41  3.13595854332114d4, 4.35061725080755d2/
42  data (a2(i),i=1,12)/-4.46620341924942d-15, -1.58654991146236d-12, &
43  -4.44467627042232d-10, -6.84738791621745d-8, &
44  -6.64932238528105d-6, -3.69976170193942d-4, &
45  -1.12295393687006d-2, -1.60926102124442d-1, &
46  -8.52408612877447d-1, -7.45519953763928d-1, &
47  2.98435207466372d0, 1.0d0/
48  data (b2(i),i=1,12)/-2.23310170962369d-15, -7.94193282071464d-13, &
49  -2.22564376956228d-10, -3.43299431079845d-8, &
50  -3.33919612678907d-6, -1.86432212187088d-4, &
51  -5.69764436880529d-3, -8.34904593067194d-2, &
52  -4.78770844009440d-1, -4.99759250374148d-1, &
53  1.86795964993052d0, 4.16485970495288d-1/
54 
55 
56  if (x .lt. 2.0d0) then
57  xx = exp(x)
58  rn = xx + a1(m1)
59  do i=m1-1,1,-1
60  rn = rn*xx + a1(i)
61  enddo
62  den = b1(k1+1)
63  do i=k1,1,-1
64  den = den*xx + b1(i)
65  enddo
66  zfermim12 = xx * rn/den
67 !
68  else
69  xx = 1.0d0/(x*x)
70  rn = xx + a2(m2)
71  do i=m2-1,1,-1
72  rn = rn*xx + a2(i)
73  enddo
74  den = b2(k2+1)
75  do i=k2,1,-1
76  den = den*xx + b2(i)
77  enddo
78  zfermim12 = sqrt(x)*rn/den
79  end if
80  return
81  end function zfermim12
82 
83 
84 
85 
86 
87 
88  real*8 function zfermi12(x)
89  include 'implno.dek'
90 !
91 ! this routine applies a rational function expansion to get the fermi
92 ! function of order 1/2 evaluated at x. maximum error is 5.47d-13.
93 ! reference: antia apjs 84,101 1993
94 !
95 ! declare
96  integer :: i,m1,k1,m2,k2
97  real*8 :: x,an,a1(12),b1(12),a2(12),b2(12),rn,den,xx
98 
99 
100 ! load the coefficients of the expansion
101  data an,m1,k1,m2,k2 /0.5d0, 7, 7, 10, 11/
102  data (a1(i),i=1,8)/5.75834152995465d6, 1.30964880355883d7, &
103  1.07608632249013d7, 3.93536421893014d6, &
104  6.42493233715640d5, 4.16031909245777d4, &
105  7.77238678539648d2, 1.0d0/
106  data (b1(i),i=1,8)/6.49759261942269d6, 1.70750501625775d7, &
107  1.69288134856160d7, 7.95192647756086d6, &
108  1.83167424554505d6, 1.95155948326832d5, &
109  8.17922106644547d3, 9.02129136642157d1/
110  data (a2(i),i=1,11)/4.85378381173415d-14, 1.64429113030738d-11, &
111  3.76794942277806d-9, 4.69233883900644d-7, &
112  3.40679845803144d-5, 1.32212995937796d-3, &
113  2.60768398973913d-2, 2.48653216266227d-1, &
114  1.08037861921488d0, 1.91247528779676d0, &
115  1.0d0/
116  data (b2(i),i=1,12)/7.28067571760518d-14, 2.45745452167585d-11, &
117  5.62152894375277d-9, 6.96888634549649d-7, &
118  5.02360015186394d-5, 1.92040136756592d-3, &
119  3.66887808002874d-2, 3.24095226486468d-1, &
120  1.16434871200131d0, 1.34981244060549d0, &
121  2.01311836975930d-1, -2.14562434782759d-2/
122 
123 
124  if (x .lt. 2.0d0) then
125  xx = exp(x)
126  rn = xx + a1(m1)
127  do i=m1-1,1,-1
128  rn = rn*xx + a1(i)
129  enddo
130  den = b1(k1+1)
131  do i=k1,1,-1
132  den = den*xx + b1(i)
133  enddo
134  zfermi12 = xx * rn/den
135 
136  else
137  xx = 1.0d0/(x*x)
138  rn = xx + a2(m2)
139  do i=m2-1,1,-1
140  rn = rn*xx + a2(i)
141  enddo
142  den = b2(k2+1)
143  do i=k2,1,-1
144  den = den*xx + b2(i)
145  enddo
146  zfermi12 = x*sqrt(x)*rn/den
147  end if
148  return
149  end function zfermi12
150 
151 
152 
153 
154 
155  real*8 function zfermi1(x)
156  include 'implno.dek'
157 !
158 ! this routine applies a rational function expansion to get the fermi
159 ! function of order 1 evaluated at x. maximum error is 1.0e-8.
160 ! reference: antia priv comm. 11sep94
161 !
162 ! declare
163  integer :: i,m1,k1,m2,k2
164  real*8 :: x,an,a1(12),b1(12),a2(12),b2(12),rn,den,xx
165 
166 ! load the coefficients of the expansion
167  data an,m1,k1,m2,k2 /1.0, 7, 4, 9, 5/
168  data (a1(i),i=1,8)/-7.606458638543d7, -1.143519707857d8, &
169  -5.167289383236d7, -7.304766495775d6, &
170  -1.630563622280d5, 3.145920924780d3, &
171  -7.156354090495d1, 1.0d0/
172  data (b1(i),i=1,5)/-7.606458639561d7, -1.333681162517d8, &
173  -7.656332234147d7, -1.638081306504d7, &
174  -1.044683266663d6/
175  data (a2(i),i=1,10)/-3.493105157219d-7, -5.628286279892d-5, &
176  -5.188757767899d-3, -2.097205947730d-1, &
177  -3.353243201574d0, -1.682094530855d1, &
178  -2.042542575231d1, 3.551366939795d0, &
179  -2.400826804233d0, 1.0d0/
180  data (b2(i),i=1,6)/-6.986210315105d-7, -1.102673536040d-4, &
181  -1.001475250797d-2, -3.864923270059d-1, &
182  -5.435619477378d0, -1.563274262745d1/
183 
184 
185  if (x .lt. 2.0d0) then
186  xx = exp(x)
187  rn = xx + a1(m1)
188  do i=m1-1,1,-1
189  rn = rn*xx + a1(i)
190  enddo
191  den = b1(k1+1)
192  do i=k1,1,-1
193  den = den*xx + b1(i)
194  enddo
195  zfermi1 = xx * rn/den
196 
197  else
198  xx = 1.0d0/(x*x)
199  rn = xx + a2(m2)
200  do i=m2-1,1,-1
201  rn = rn*xx + a2(i)
202  enddo
203  den = b2(k2+1)
204  do i=k2,1,-1
205  den = den*xx + b2(i)
206  enddo
207  zfermi1 = x*x*rn/den
208  end if
209  return
210  end function zfermi1
211 
212 
213 
214 
215 
216  real*8 function zfermi32(x)
217  include 'implno.dek'
218 !
219 ! this routine applies a rational function expansion to get the fermi
220 ! function of order 3/2 evaluated at x. maximum error is 5.07d-13.
221 ! reference: antia apjs 84,101 1993
222 !
223 ! declare
224  integer :: i,m1,k1,m2,k2
225  real*8 :: x,an,a1(12),b1(12),a2(12),b2(12),rn,den,xx
226 
227 ! load the coefficients of the expansion
228  data an,m1,k1,m2,k2 /1.5d0, 6, 7, 9, 10/
229  data (a1(i),i=1,7)/4.32326386604283d4, 8.55472308218786d4, &
230  5.95275291210962d4, 1.77294861572005d4, &
231  2.21876607796460d3, 9.90562948053193d1, &
232  1.0d0/
233  data (b1(i),i=1,8)/3.25218725353467d4, 7.01022511904373d4, &
234  5.50859144223638d4, 1.95942074576400d4, &
235  3.20803912586318d3, 2.20853967067789d2, &
236  5.05580641737527d0, 1.99507945223266d-2/
237  data (a2(i),i=1,10)/2.80452693148553d-13, 8.60096863656367d-11, &
238  1.62974620742993d-8, 1.63598843752050d-6, &
239  9.12915407846722d-5, 2.62988766922117d-3, &
240  3.85682997219346d-2, 2.78383256609605d-1, &
241  9.02250179334496d-1, 1.0d0/
242  data (b2(i),i=1,11)/7.01131732871184d-13, 2.10699282897576d-10, &
243  3.94452010378723d-8, 3.84703231868724d-6, &
244  2.04569943213216d-4, 5.31999109566385d-3, &
245  6.39899717779153d-2, 3.14236143831882d-1, &
246  4.70252591891375d-1, -2.15540156936373d-2, &
247  2.34829436438087d-3/
248 
249 
250  if (x .lt. 2.0d0) then
251  xx = exp(x)
252  rn = xx + a1(m1)
253  do i=m1-1,1,-1
254  rn = rn*xx + a1(i)
255  enddo
256  den = b1(k1+1)
257  do i=k1,1,-1
258  den = den*xx + b1(i)
259  enddo
260  zfermi32 = xx * rn/den
261 
262  else
263  xx = 1.0d0/(x*x)
264  rn = xx + a2(m2)
265  do i=m2-1,1,-1
266  rn = rn*xx + a2(i)
267  enddo
268  den = b2(k2+1)
269  do i=k2,1,-1
270  den = den*xx + b2(i)
271  enddo
272  zfermi32 = x*x*sqrt(x)*rn/den
273  end if
274  return
275  end function zfermi32
276 
277 
278 
279 
280  real*8 function zfermi2(x)
281  include 'implno.dek'
282 !
283 ! this routine applies a rational function expansion to get the fermi
284 ! function of order 2 evaluated at x. maximum error is 1.0e-8.
285 ! reference: antia priv comm. 11sep94
286 !
287 ! declare
288  integer :: i,m1,k1,m2,k2
289  real*8 :: x,an,a1(12),b1(12),a2(12),b2(12),rn,den,xx
290 
291 ! load the coefficients of the expansion
292  data an,m1,k1,m2,k2 /2.0, 7, 4, 5, 9/
293  data (a1(i),i=1,8)/-1.434885992395d8, -2.001711155617d8, &
294  -8.507067153428d7, -1.175118281976d7, &
295  -3.145120854293d5, 4.275771034579d3, &
296  -8.069902926891d1, 1.0d0/
297  data (b1(i),i=1,5)/-7.174429962316d7, -1.090535948744d8, &
298  -5.350984486022d7, -9.646265123816d6, &
299  -5.113415562845d5/
300  data (a2(i),i=1,6)/ 6.919705180051d-8, 1.134026972699d-5, &
301  7.967092675369d-4, 2.432500578301d-2, &
302  2.784751844942d-1, 1.0d0/
303  data (b2(i),i=1,10)/ 2.075911553728d-7, 3.197196691324d-5, &
304  2.074576609543d-3, 5.250009686722d-2, &
305  3.171705130118d-1, -1.147237720706d-1, &
306  6.638430718056d-2, -1.356814647640d-2, &
307  -3.648576227388d-2, 3.621098757460d-2/
308 
309 
310  if (x .lt. 2.0d0) then
311  xx = exp(x)
312  rn = xx + a1(m1)
313  do i=m1-1,1,-1
314  rn = rn*xx + a1(i)
315  enddo
316  den = b1(k1+1)
317  do i=k1,1,-1
318  den = den*xx + b1(i)
319  enddo
320  zfermi2 = xx * rn/den
321 
322  else
323  xx = 1.0d0/(x*x)
324  rn = xx + a2(m2)
325  do i=m2-1,1,-1
326  rn = rn*xx + a2(i)
327  enddo
328  den = b2(k2+1)
329  do i=k2,1,-1
330  den = den*xx + b2(i)
331  enddo
332  zfermi2 = x*x*x*rn/den
333  end if
334  return
335  end function zfermi2
336 
337 
338 
339 
340 
341 
342  real*8 function zfermi52(x)
343  include 'implno.dek'
344 !
345 ! this routine applies a rational function expansion to get the fermi
346 ! function of order 5/2 evaluated at x. maximum error is 2.47d-13.
347 ! reference: antia apjs 84,101 1993
348 !
349 ! declare
350  integer :: i,m1,k1,m2,k2
351  real*8 :: x,an,a1(12),b1(12),a2(12),b2(12),rn,den,xx
352 
353 ! load the coefficients of the expansion
354  data an,m1,k1,m2,k2 /2.5d0, 6, 7, 10, 9/
355  data (a1(i),i=1,7)/6.61606300631656d4, 1.20132462801652d5, &
356  7.67255995316812d4, 2.10427138842443d4, &
357  2.44325236813275d3, 1.02589947781696d2, &
358  1.0d0/
359  data (b1(i),i=1,8)/1.99078071053871d4, 3.79076097261066d4, &
360  2.60117136841197d4, 7.97584657659364d3, &
361  1.10886130159658d3, 6.35483623268093d1, &
362  1.16951072617142d0, 3.31482978240026d-3/
363  data (a2(i),i=1,11)/8.42667076131315d-12, 2.31618876821567d-9, &
364  3.54323824923987d-7, 2.77981736000034d-5, &
365  1.14008027400645d-3, 2.32779790773633d-2, &
366  2.39564845938301d-1, 1.24415366126179d0, &
367  3.18831203950106d0, 3.42040216997894d0, &
368  1.0d0/
369  data (b2(i),i=1,10)/2.94933476646033d-11, 7.68215783076936d-9, &
370  1.12919616415947d-6, 8.09451165406274d-5, &
371  2.81111224925648d-3, 3.99937801931919d-2, &
372  2.27132567866839d-1, 5.31886045222680d-1, &
373  3.70866321410385d-1, 2.27326643192516d-2/
374 
375 
376  if (x .lt. 2.0d0) then
377  xx = exp(x)
378  rn = xx + a1(m1)
379  do i=m1-1,1,-1
380  rn = rn*xx + a1(i)
381  enddo
382  den = b1(k1+1)
383  do i=k1,1,-1
384  den = den*xx + b1(i)
385  enddo
386  zfermi52 = xx * rn/den
387 
388  else
389  xx = 1.0d0/(x*x)
390  rn = xx + a2(m2)
391  do i=m2-1,1,-1
392  rn = rn*xx + a2(i)
393  enddo
394  den = b2(k2+1)
395  do i=k2,1,-1
396  den = den*xx + b2(i)
397  enddo
398  zfermi52 = x*x*x*sqrt(x)*rn/den
399  end if
400  return
401  end function zfermi52
402 
403 
404 
405 
406 
407  real*8 function zfermi3(x)
408  include 'implno.dek'
409 !
410 ! this routine applies a rational function expansion to get the fermi
411 ! function of order 3 evaluated at x. maximum error is 1.0e-8.
412 ! reference: antia priv comm. 11sep94
413 !
414 ! declare
415  integer :: i,m1,k1,m2,k2
416  real*8 :: x,an,a1(12),b1(12),a2(12),b2(12),rn,den,xx
417 
418 ! load the coefficients of the expansion
419  data an,m1,k1,m2,k2 /3.0, 4, 6, 7, 7/
420  data (a1(i),i=1,5)/ 6.317036716422d2, 7.514163924637d2, &
421  2.711961035750d2, 3.274540902317d1, &
422  1.0d0/
423  data (b1(i),i=1,7)/ 1.052839452797d2, 1.318163114785d2, &
424  5.213807524405d1, 7.500064111991d0, &
425  3.383020205492d-1, 2.342176749453d-3, &
426  -8.445226098359d-6/
427  data (a2(i),i=1,8)/ 1.360999428425d-8, 1.651419468084d-6, &
428  1.021455604288d-4, 3.041270709839d-3, &
429  4.584298418374d-2, 3.440523212512d-1, &
430  1.077505444383d0, 1.0d0/
431  data (b2(i),i=1,8)/ 5.443997714076d-8, 5.531075760054d-6, &
432  2.969285281294d-4, 6.052488134435d-3, &
433  5.041144894964d-2, 1.048282487684d-1, &
434  1.280969214096d-2, -2.851555446444d-3/
435 
436 
437  if (x .lt. 2.0d0) then
438  xx = exp(x)
439  rn = xx + a1(m1)
440  do i=m1-1,1,-1
441  rn = rn*xx + a1(i)
442  enddo
443  den = b1(k1+1)
444  do i=k1,1,-1
445  den = den*xx + b1(i)
446  enddo
447  zfermi3 = xx * rn/den
448 
449  else
450  xx = 1.0d0/(x*x)
451  rn = xx + a2(m2)
452  do i=m2-1,1,-1
453  rn = rn*xx + a2(i)
454  enddo
455  den = b2(k2+1)
456  do i=k2,1,-1
457  den = den*xx + b2(i)
458  enddo
459  zfermi3 = x*x*x*x*rn/den
460  end if
461  return
462  end function zfermi3
463 
464 
465 
466 
467 
468  real*8 function ifermim12(f)
469  include 'implno.dek'
470 !
471 ! this routine applies a rational function expansion to get the inverse
472 ! fermi function of order -1/2 when it is equal to f.
473 ! maximum error is 3.03d-9. reference: antia apjs 84,101 1993
474 !
475 ! declare
476  integer :: i,m1,k1,m2,k2
477  real*8 :: f,an,a1(12),b1(12),a2(12),b2(12),rn,den,ff
478 
479 ! load the coefficients of the expansion
480  data an,m1,k1,m2,k2 /-0.5d0, 5, 6, 6, 6/
481  data (a1(i),i=1,6)/-1.570044577033d4, 1.001958278442d4, &
482  -2.805343454951d3, 4.121170498099d2, &
483  -3.174780572961d1, 1.0d0/
484  data (b1(i),i=1,7)/-2.782831558471d4, 2.886114034012d4, &
485  -1.274243093149d4, 3.063252215963d3, &
486  -4.225615045074d2, 3.168918168284d1, &
487  -1.008561571363d0/
488  data (a2(i),i=1,7)/ 2.206779160034d-8, -1.437701234283d-6, &
489  6.103116850636d-5, -1.169411057416d-3, &
490  1.814141021608d-2, -9.588603457639d-2, &
491  1.0d0/
492  data (b2(i),i=1,7)/ 8.827116613576d-8, -5.750804196059d-6, &
493  2.429627688357d-4, -4.601959491394d-3, &
494  6.932122275919d-2, -3.217372489776d-1, &
495  3.124344749296d0/
496 
497  if (f .lt. 4.0d0) then
498  rn = f + a1(m1)
499  do i=m1-1,1,-1
500  rn = rn*f + a1(i)
501  enddo
502  den = b1(k1+1)
503  do i=k1,1,-1
504  den = den*f + b1(i)
505  enddo
506  ifermim12 = log(f * rn/den)
507 
508  else
509  ff = 1.0d0/f**(1.0d0/(1.0d0 + an))
510  rn = ff + a2(m2)
511  do i=m2-1,1,-1
512  rn = rn*ff + a2(i)
513  enddo
514  den = b2(k2+1)
515  do i=k2,1,-1
516  den = den*ff + b2(i)
517  enddo
518  ifermim12 = rn/(den*ff)
519  end if
520  return
521  end function ifermim12
522 
523 
524 
525 
526 
527 
528  real*8 function ifermi12(f)
529  include 'implno.dek'
530 !
531 ! this routine applies a rational function expansion to get the inverse
532 ! fermi function of order 1/2 when it is equal to f.
533 ! maximum error is 4.19d-9. reference: antia apjs 84,101 1993
534 !
535 ! declare
536  integer :: i,m1,k1,m2,k2
537  real*8 :: f,an,a1(12),b1(12),a2(12),b2(12),rn,den,ff
538 
539 ! load the coefficients of the expansion
540  data an,m1,k1,m2,k2 /0.5d0, 4, 3, 6, 5/
541  data (a1(i),i=1,5)/ 1.999266880833d4, 5.702479099336d3, &
542  6.610132843877d2, 3.818838129486d1, &
543  1.0d0/
544  data (b1(i),i=1,4)/ 1.771804140488d4, -2.014785161019d3, &
545  9.130355392717d1, -1.670718177489d0/
546  data (a2(i),i=1,7)/-1.277060388085d-2, 7.187946804945d-2, &
547  -4.262314235106d-1, 4.997559426872d-1, &
548  -1.285579118012d0, -3.930805454272d-1, &
549  1.0d0/
550  data (b2(i),i=1,6)/-9.745794806288d-3, 5.485432756838d-2, &
551  -3.299466243260d-1, 4.077841975923d-1, &
552  -1.145531476975d0, -6.067091689181d-2/
553 
554 
555  if (f .lt. 4.0d0) then
556  rn = f + a1(m1)
557  do i=m1-1,1,-1
558  rn = rn*f + a1(i)
559  enddo
560  den = b1(k1+1)
561  do i=k1,1,-1
562  den = den*f + b1(i)
563  enddo
564  ifermi12 = log(f * rn/den)
565 
566  else
567  ff = 1.0d0/f**(1.0d0/(1.0d0 + an))
568  rn = ff + a2(m2)
569  do i=m2-1,1,-1
570  rn = rn*ff + a2(i)
571  enddo
572  den = b2(k2+1)
573  do i=k2,1,-1
574  den = den*ff + b2(i)
575  enddo
576  ifermi12 = rn/(den*ff)
577  end if
578  return
579  end function ifermi12
580 
581 
582 
583 
584 
585 
586  real*8 function ifermi32(f)
587  include 'implno.dek'
588 !
589 ! this routine applies a rational function expansion to get the inverse
590 ! fermi function of order 3/2 when it is equal to f.
591 ! maximum error is 2.26d-9. reference: antia apjs 84,101 1993
592 !
593 ! declare
594  integer :: i,m1,k1,m2,k2
595  real*8 :: f,an,a1(12),b1(12),a2(12),b2(12),rn,den,ff
596 
597 ! load the coefficients of the expansion
598  data an,m1,k1,m2,k2 /1.5d0, 3, 4, 6, 5/
599  data (a1(i),i=1,4)/ 1.715627994191d2, 1.125926232897d2, &
600  2.056296753055d1, 1.0d0/
601  data (b1(i),i=1,5)/ 2.280653583157d2, 1.193456203021d2, &
602  1.167743113540d1, -3.226808804038d-1, &
603  3.519268762788d-3/
604  data (a2(i),i=1,7)/-6.321828169799d-3, -2.183147266896d-2, &
605  -1.057562799320d-1, -4.657944387545d-1, &
606  -5.951932864088d-1, 3.684471177100d-1, &
607  1.0d0/
608  data (b2(i),i=1,6)/-4.381942605018d-3, -1.513236504100d-2, &
609  -7.850001283886d-2, -3.407561772612d-1, &
610  -5.074812565486d-1, -1.387107009074d-1/
611 
612 
613  if (f .lt. 4.0d0) then
614  rn = f + a1(m1)
615  do i=m1-1,1,-1
616  rn = rn*f + a1(i)
617  enddo
618  den = b1(k1+1)
619  do i=k1,1,-1
620  den = den*f + b1(i)
621  enddo
622  ifermi32 = log(f * rn/den)
623 
624  else
625  ff = 1.0d0/f**(1.0d0/(1.0d0 + an))
626  rn = ff + a2(m2)
627  do i=m2-1,1,-1
628  rn = rn*ff + a2(i)
629  enddo
630  den = b2(k2+1)
631  do i=k2,1,-1
632  den = den*ff + b2(i)
633  enddo
634  ifermi32 = rn/(den*ff)
635  end if
636  return
637  end function ifermi32
638 
639 
640 
641 
642 
643  real*8 function ifermi52(f)
644  include 'implno.dek'
645 !
646 ! this routine applies a rational function expansion to get the inverse
647 ! fermi function of order 5/2 when it is equal to f.
648 ! maximum error is 6.17d-9. reference: antia apjs 84,101 1993
649 !
650 ! declare
651  integer :: i,m1,k1,m2,k2
652  real*8 :: f,an,a1(12),b1(12),a2(12),b2(12),rn,den,ff
653 
654 ! load the coefficients of the expansion
655  data an,m1,k1,m2,k2 /2.5d0, 2, 3, 6, 6/
656  data (a1(i),i=1,3)/ 2.138969250409d2, 3.539903493971d1, &
657  1.0d0/
658  data (b1(i),i=1,4)/ 7.108545512710d2, 9.873746988121d1, &
659  1.067755522895d0, -1.182798726503d-2/
660  data (a2(i),i=1,7)/-3.312041011227d-2, 1.315763372315d-1, &
661  -4.820942898296d-1, 5.099038074944d-1, &
662  5.495613498630d-1, -1.498867562255d0, &
663  1.0d0/
664  data (b2(i),i=1,7)/-2.315515517515d-2, 9.198776585252d-2, &
665  -3.835879295548d-1, 5.415026856351d-1, &
666  -3.847241692193d-1, 3.739781456585d-2, &
667  -3.008504449098d-2/
668 
669 
670  if (f .lt. 4.0d0) then
671  rn = f + a1(m1)
672  do i=m1-1,1,-1
673  rn = rn*f + a1(i)
674  enddo
675  den = b1(k1+1)
676  do i=k1,1,-1
677  den = den*f + b1(i)
678  enddo
679  ifermi52 = log(f * rn/den)
680 
681  else
682  ff = 1.0d0/f**(1.0d0/(1.0d0 + an))
683  rn = ff + a2(m2)
684  do i=m2-1,1,-1
685  rn = rn*ff + a2(i)
686  enddo
687  den = b2(k2+1)
688  do i=k2,1,-1
689  den = den*ff + b2(i)
690  enddo
691  ifermi52 = rn/(den*ff)
692  end if
693  return
694  end function ifermi52
zfermi1
real *8 function zfermi1(x)
Definition: funct_fermi1.f90:156
zfermi3
real *8 function zfermi3(x)
Definition: funct_fermi1.f90:408
ifermi32
real *8 function ifermi32(f)
Definition: funct_fermi1.f90:587
ifermi12
real *8 function ifermi12(f)
Definition: funct_fermi1.f90:529
ifermi52
real *8 function ifermi52(f)
Definition: funct_fermi1.f90:644
zfermi32
real *8 function zfermi32(x)
Definition: funct_fermi1.f90:217
zfermi12
real *8 function zfermi12(x)
Definition: funct_fermi1.f90:89
zfermi52
real *8 function zfermi52(x)
Definition: funct_fermi1.f90:343
zfermim12
real *8 function zfermim12(x)
Definition: funct_fermi1.f90:22
zfermi2
real *8 function zfermi2(x)
Definition: funct_fermi1.f90:281
ifermim12
real *8 function ifermim12(f)
Definition: funct_fermi1.f90:469