File: snell.f

package info (click to toggle)
scilab 4.0-12
  • links: PTS
  • area: non-free
  • in suites: etch, etch-m68k
  • size: 100,640 kB
  • ctags: 57,333
  • sloc: ansic: 377,889; fortran: 242,862; xml: 179,819; tcl: 42,062; sh: 10,593; ml: 9,441; makefile: 4,377; cpp: 1,354; java: 621; csh: 260; yacc: 247; perl: 130; lex: 126; asm: 72; lisp: 30
file content (53 lines) | stat: -rw-r--r-- 1,040 bytes parent folder | download | duplicates (14) 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
 
C/MEMBR ADD NAME=SNELL,SSI=0
      subroutine snell(dsn2,du, dk, dq)
c!purpose
c calculation of the jacobi's elliptic function sn(u,k)
c
c external calculation of the parameter necessary
c dk = k($k)
c dq = exp(-pi*k'/k) ... (jacobi's nome)
c!
c
      double precision dpi, domi
      double precision de, dz, dpi2, dq, dm, du, dk, dc, dqq, dh, dq1,
     *    dq2
c
      domi=2.0d+0*dlamch('p')
      dpi=4.0d+0*atan(1.0d+0)
c
      data de, dz /1.0d+0,2.0d+0/
c
      dpi2 = dpi/dz
      if (abs(dq).ge.de) go to 30
c
      dm = dpi2*du/dk
      dc = dz*dm
      dc =  cos(dc)
c
      dm =  sin(dm)*dk/dpi2
      dqq = dq*dq
      dq1 = dq
      dq2 = dqq
c
      do 10 i=1,100
        dh = (de-dq1)/(de-dq2)
        dh = dh*dh
        dh = dh*(de-dz*dq2*dc+dq2*dq2)
        dh = dh/(de-dz*dq1*dc+dq1*dq1)
        dm = dm*dh
c
        dh = abs(de-dh)
        if (dh.lt.domi) go to 20
c
        dq1 = dq1*dqq
        dq2 = dq2*dqq
  10  continue
c
      go to 30
c
  20  dsn2 = dm
      return
c
  30  dsn2 = 0.0d+0
      return
      end