V3FIT
surfaces_plot.f
1  subroutine surfaces_plot(g1,g2,device)
2  use precision
3  use vname0
4  use gfile
5  use mapout
6  implicit none
7  type(geqdsk) :: g1
8  type(rszs) :: g2
9  character*(*) :: device
10  real*4, dimension(:), allocatable :: sqphin,thetap,xp,yp
11  real, pointer :: xpnt(:,:),ypnt(:,:)
12  real(rprec) :: pi
13  real*4 xmin, xmax, zmin, zmax
14  integer :: id0, id1, id, pgopen
15  integer :: jpsi, itht, j, ij, m
16  character(len=60) :: mlabel
17  interface
18  subroutine evals_boundary(g2,xp,yp,icount)
19  use vname0
20  use mapout
21  real(rprec), dimension(0:mu-1,0:0) :: rbc,zbs,rbs,zbc
22  real*4, dimension(:), allocatable :: rb, zb
23  real*4, dimension(*) :: xp, yp
24  type(rszs) :: g2
25  end subroutine evals_boundary
26  end interface
27  pi = 4 * atan(1._dbl)
28  call pgask(.true.)
29  zmin=-165. ! allow for D3D coils
30  zmax=165.
31  xmin=80.
32  xmax=270.
33 ! call pgpap (11., (zmax-zmin)/(xmax-xmin)/2) !
34 ! if(INDEX(TRIM(device),'p')/=0)
35  call pgpap (8., (zmax-zmin)/(xmax-xmin)/2)
36  call pgsubp(2,1)
37  call pgsch (2.0) ! set height of text
38  call pgsci(1.)
39  call pgenv (xmin, xmax, zmin, zmax,1,0) ! define the subplot area
40  call pglab('R(cm)','Z(cm)','mapcode')
41  call pgscf (1 ) ! select simple typeface
42  call pgslw(5)
43  m=g1%limitr;allocate(xp(m),yp(m))
44  xp(1:m)=100*g1%xlim(1:m);yp(1:m)=100*g1%ylim(1:m)
45  call pgline(m,xp,yp)
46  if( allocated(xp) )deallocate(xp,yp)
47  allocate(xp(nu),yp(nu))
48  xp=0;yp=0
49  call evals_boundary(g2,xp,yp,m)
50  xp=xp*100;yp=yp*100
51  call pgpt(m,xp,yp,5)
52  if( allocated(xp) )deallocate(xp,yp)
53  call pltcol
54  call pgslw(2)
55  call pgsci(3)
56  m=g1%nbdry;allocate(xp(m),yp(m))
57  xp(1:m)=100*g1%rbdry(1:m);yp(1:m)=100*g1%zbdry(1:m)
58  call pgpt(m,xp,yp,3)
59  if( allocated(xp) )deallocate(xp,yp)
60  call pgslw(2)
61  itht=g2%nthet
62  jpsi=g2%npsi
63  allocate(xp(g2%nthet),yp(g2%nthet))
64  call pgsci(4)
65  do j=1,jpsi,20
66  xp(1:itht)=100*g2%rs(j,1:itht)
67  yp(1:itht)=100*g2%zs(j,1:itht)
68  call pgline(itht,xp,yp)
69  enddo
70  if( allocated(xp) )deallocate(xp,yp)
71  allocate(xp(g2%npsi),yp(g2%npsi))
72  call pgsci(5)
73  ij=itht/16
74  do j=1,itht,ij
75  xp(1:jpsi)=100*g2%rs(1:jpsi,j)
76  yp(1:jpsi)=100*g2%zs(1:jpsi,j)
77  call pgline(jpsi,xp,yp)
78  enddo
79  call pgsci(2)
80  j=1
81  xp(1:jpsi)=100*g2%rs(1:jpsi,j)
82  yp(1:jpsi)=100*g2%zs(1:jpsi,j)
83  call pgline(jpsi,xp,yp)
84  j=itht/2
85  xp(1:jpsi)=100*g2%rs(1:jpsi,j)
86  yp(1:jpsi)=100*g2%zs(1:jpsi,j)
87  call pgline(jpsi,xp,yp)
88  call pgsci(1)
89 c second plot
90  call pgenv (xmin, xmax, zmin, zmax,1,0) ! define the subplot area
91  call pglab('R(cm)','Z(cm)','DESCUR v mapped eqdsk')
92  write(mlabel,fmt='("% boundary flux=",f8.4)')
93  . 100.*g2%fraction_bndry
94  CALL pgmtxt('B',2.25,0.5,0.5,trim(mlabel))
95  write(mlabel,fmt='("mpol=",i3.3)')mu
96  CALL pgmtxt('T',0.325,0.5,0.5,trim(mlabel))
97  call pgscf (1 ) ! select simple typeface
98  call pgsch (2.0) ! set height of text
99  call pgslw(5)
100  if( allocated(xp) )deallocate(xp,yp)
101  m=g1%limitr;allocate(xp(m),yp(m))
102  xp(1:m)=100*g1%xlim(1:m);yp(1:m)=100*g1%ylim(1:m)
103  call pgline(m,xp,yp)
104  if( allocated(xp) )deallocate(xp,yp)
105  call pltcol
106  call pgslw(2)
107  if( allocated(xp) )deallocate(xp,yp)
108  allocate(xp(g2%nthet),yp(g2%nthet))
109  call pgsci(4)
110  j=g2%npsi
111  xp(1:itht)=100*g2%rs(j,1:itht)
112  yp(1:itht)=100*g2%zs(j,1:itht)
113  call pgpt(itht,xp,yp,5)
114  if( allocated(xp) )deallocate(xp,yp)
115  call pgsci(3)
116  if( allocated(xp) )deallocate(xp,yp)
117  allocate(xp(nu),yp(nu))
118  xp=0;yp=0
119  call evals_boundary(g2,xp,yp,m)
120  xp=xp*100;yp=yp*100
121  call pgslw(5)
122  call pgsci(2)
123  call pgline(m,xp,yp)
124  end subroutine surfaces_plot
125  subroutine pltcol
126  parameter(ncoil=18)
127  dimension rc(ncoil), zc(ncoil), wc(ncoil),
128  .hc(ncoil),xx(5),yy(5)
129  dimension ac(ncoil),ac2(ncoil)
130 c
131  data rc/
132  &.8608,.8614,.8628,.8611,1.0041,2.6124,
133  &2.3733,1.2518,1.6890,.8608,.8607,.8611,
134  &.8630,1.0025,2.6124,2.3834,1.2524,1.6889/
135  data zc/
136  &.16830,.50810,.84910,1.1899,1.5169,0.4376,
137  &1.1171,1.6019,1.5874,-.17370,-.51350,-.85430,
138  &-1.1957,-1.5169,-0.4376,-1.1171,-1.6027,-1.5780/
139  data wc/
140  &.0508,.0508,.0508,.0508,.13920,0.17320,
141  &0.1880,.23490,.16940,.0508,.0508,.0508,
142  &.0508,.13920,0.17320,0.1880,.23490,.16940/
143  data hc/
144  &.32106,.32106,.32106,.32106,.11940,0.1946,
145  &0.16920,.08510,.13310,.32106,.32106,.32106,
146  &.32106,.11940,0.1946,0.16920,.08510,.13310/
147  data ac/
148  &0.,0.,0.,0.,45.0,0.,
149  &0.,0.,0.,0.,0.,0.,
150  &0.,-45.0,0.,0.,0.,0./
151  data ac2/
152  &0.,0.,0.,0.,0.,92.40,
153  &108.06,0.,0.,0.,0.,0.,
154  &0.,0.,-92.40,-108.06,0.,0./
155  do i=1,ncoil
156  x= rc(i)
157  y= zc(i)
158  dx=wc(i)/2.
159  dy=hc(i)/2.
160 c gfortran does not provide the sind and cosd functions.
161  sn1=sin(ac(i)*pi/180)
162  cos1=cos(ac(i)*pi/180)
163  tac=sn1/cos1
164  sn2=sin(ac2(i)*pi/180)
165  cos2=cos(ac2(i)*pi/180)
166  tac2=sn2/cos2
167  if(ac2(i).eq.0) go to 40
168  xx(1)=x-dx-dy/tac2
169  xx(3)=x+dx+dy/tac2
170  xx(5)=xx(1)
171  xx(2)=xx(3)-wc(i)
172  xx(4)=xx(1)+wc(i)
173 c
174  yy(1)=y-dy
175  yy(2)=y+dy
176  yy(3)=y+dy
177  yy(4)=y-dy
178  yy(5)=y-dy
179  go to 50
180  40 continue
181 c
182  dac=0.5*wc(i)*tac
183  xx(1)=x-dx
184  xx(2)=x-dx
185  xx(3)=x+dx
186  xx(4)=x+dx
187  xx(5)=x-dx
188 c
189  yy(1)=y-dy-dac
190  yy(2)=y+dy-dac
191  yy(3)=y+dy+dac
192  yy(4)=y-dy+dac
193  yy(5)=y-dy-dac
194  50 continue
195  do is=1,5
196  xx(is)=xx(is)*100
197  yy(is)=yy(is)*100
198  enddo
199  call pgsci(1)
200  call pgslw(3)
201  call pgline(5,xx,yy)
202  call pgsci(2)
203  call pgpoly(5,xx,yy)
204  enddo
205  end subroutine pltcol
206  subroutine arc_integral(rv,zv,jpsi,itht,arcl)
207  real arcl(jpsi-1)
208  real, dimension(jpsi,itht) :: rv,zv
209  real, dimension(:), allocatable :: dl
210  if(.not.allocated(dl))allocate(dl(itht-1))
211  pi=4 * atan(1.0)
212  itht=itht
213  do k=1,jpsi
214  dl=sqrt( (rv(k,2:itht)-rv(k,1:itht-1))**2 +
215  & (zv(k,2:itht)-zv(k,1:itht-1))**2 )
216  arcl(k)=sum(dl)
217  enddo
218  deallocate(dl)
219  end subroutine arc_integral