Line | |
---|
1 | subroutine quad(zn) |
---|
2 | c--------------------------------------------------------------------------- |
---|
3 | save |
---|
4 | c |
---|
5 | include 'param_sz.h' |
---|
6 | include 'constcom.h' |
---|
7 | include 'pcordcom.h' |
---|
8 | include 'syscom.h' |
---|
9 | include 'ucom.h' |
---|
10 | c |
---|
11 | dimension e(6) |
---|
12 | c-------------------------------------------------------------------------- |
---|
13 | c* |
---|
14 | ne=rne |
---|
15 | sl=zn-z |
---|
16 | k=1 |
---|
17 | l=4 |
---|
18 | xp=bgx/bgz |
---|
19 | yp=bgy/bgz |
---|
20 | con=1.+xp**2+yp**2 |
---|
21 | z=zn |
---|
22 | caysq=el(4,ne)/(brhof*bgz) |
---|
23 | cay=sqrt(abs(caysq)) |
---|
24 | if(caysq)10,30,20 |
---|
25 | 10 continue |
---|
26 | k=4 |
---|
27 | l=1 |
---|
28 | 20 continue |
---|
29 | cayl=cay*sl |
---|
30 | s=sin(cayl) |
---|
31 | c=cos(cayl) |
---|
32 | sh=sinh(cayl) |
---|
33 | ch=cosh(cayl) |
---|
34 | e(k)=c |
---|
35 | e(k+1)=s/cay |
---|
36 | e(k+2)=-cay*s |
---|
37 | e(l)=ch |
---|
38 | e(l+1)=sh/cay |
---|
39 | e(l+2)=cay*sh |
---|
40 | xs=x |
---|
41 | ys=y |
---|
42 | x=e(1)*xs+e(2)*xp |
---|
43 | xp=e(3)*xs+e(1)*xp |
---|
44 | y=e(4)*ys+e(5)*yp |
---|
45 | yp=e(6)*ys+e(4)*yp |
---|
46 | bgx=xp*bgz |
---|
47 | bgy=yp*bgz |
---|
48 | bgz=bgz*sqrt(con/(1.+xp**2+yp**2)) |
---|
49 | return |
---|
50 | 30 continue |
---|
51 | x=x+xp*sl |
---|
52 | y=y+yp*sl |
---|
53 | return |
---|
54 | end |
---|
55 | c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* |
---|
Note: See
TracBrowser
for help on using the repository browser.