module plots ! ! (c) Copyright 1991 to 1998 by Michael A. Beer, William D. Dorland, ! P. B. Snyder, Q. P. Liu, and Gregory W. Hammett. ALL RIGHTS RESERVED. ! use itg_data use gryffin_grid use graphics use diagnostics_arrays, only : date, nplot2, mplot2 use constants implicit none integer, dimension(:), allocatable :: dash ! dash pattern to identify mode # in plots contains subroutine hmplot(a,title,tim,ipage) ! ! Plot some coefficients vs. ballooning angle for various m's. ! where a is a set of real (not complex) coefficients. ! real, dimension(:,:,:) :: a character(*) :: title real :: tim, ymin, ymax, xmin, xmax integer :: l, m, n, isympa, ipage character*32 strtime ymin=0. ymax=0. n=nplot2 do m=1,md do l=1,ldb if (a(l,m,n) > ymax) ymax=a(l,m,n) if (a(l,m,n) < ymin) ymin=a(l,m,n) enddo enddo xmin=r(1) xmax=r(ld) call setzer call plchlq(.5,.08,'ballooning coordinate',15.,0.,0.) call plchlq(.5,.97,title,20.,0.,0.) write(strtime,'(a,i3,a,f7.2)') 'n = ',nplot2,', all m''s, time = ',tim call plchlq(.98,.93,strtime,15.,0.,1.) call mapsm(xmin,xmax,ymin,ymax) n=nplot2 do m=1,md call dashdb(dash(m)) call curved(r,a(:,m,n),ldb) enddo call finish(ipage, date) return end subroutine hmplot subroutine hmplot2(field, title, tim, ipage, mplot, nplot) ! ! plots the real and imaginary mode structure of a field along the ! the magnetic field, for a single mode on a single page. ! title is a character string title for the plot ! ! Then plot the Fourier Transform of the field along B. ! ! arguments: complex, dimension(:,:,:) :: field character(*) title real tim integer mplot,nplot,ipage ! local variables: real, dimension(size(field, 1)) :: yreal, yimag, ucft, usft, work, rf complex phase,ylocal character*14 strtime character*39 strmn ! declare other local variables: real ymin,ymax,xmin,xmax integer isymspa,l,ifail interface subroutine mnmx(n,vector,vmin,vmax) integer :: n real, dimension(:) :: vector real :: vmin, vmax end subroutine mnmx end interface ! rotate the complex mode to make it purely real at r=0: phase=field(ld/2+1,mplot,nplot) if(cabs(phase) > 1.e-30)then phase=phase/cabs(phase) else phase=(1.0,0.0) endif do l=1,ld ylocal=field(l,mplot,nplot)/phase yreal(l)=real(ylocal) yimag(l)=aimag(ylocal) ! yimag(l)=-aimag(ylocal) ! We used to have to take the complex conjugate to compare with Tang ! (PF v. 23, p.2459, fig. 5), but now I realize that was because our ! sin part was the negative of the imaginary part. enddo ! debug: try a gaussian: ! do 12 l=1,ld ! yimag(l)=0.0 ! yreal(l)=exp(-(r(l)-r(ld)/2)**2/2) !12 continue ymin=0.0 ymax=0.0 call mnmx(ld,yreal,ymin,ymax) call mnmx(ld,yimag,ymin,ymax) xmin=r(1) xmax=r(ld) ! start plot call setzer ! init itext to hold up to 5 80-character lines for the legend: ! idum=linest(itext,itextsiz,80) call plchlq(.5,.08,'ballooning coordinate',15.,0.,0.) call plchlq(.5,.97,title,20.,0.,0.) write(strmn,'(a,i2,a,i2)') 'Real and Imag part vs. x, for m=', & mrr(mplot,nplot),' n=',nrr(nplot) call plchlq(.5,.93,strmn,15.,0.,0.) write(strtime,'(a,f7.2)') 'time = ',tim call plchlq(.98,.98,strtime,15.,0.,1.) call mapsm(xmin,xmax,ymin,ymax) isymspa=(ld-1)/16 call dashdb(65535) call curve(r,yreal,ld) call dashdb(21845) call curved(r,yimag,ld) call set(0.1,0.9,0.15,0.9,0.,1.,0.,1.,1) call line(.15,.1,.18,.1) call plchlq(.19,.1,'real',12.,0.,-1.) call lined(.15,.05,.18,.05) call plchlq(.19,.05,'imaginary',12.,0.,-1.) call finish(ipage, date) !*************************************************************** ! !c FFT phi(k) for plots: ! if ((qsf >= 0.0).and.(mrr(mplot,nplot) /= 0)) then ! ! F.T. of potential ! do l=1,ld-1 ! define the k_parallel grid: rf(l)=float(l-ldb/2)/(2*x0) ! 2.*r(l)/(3.14159*(ld-1))/(mrr(mplot,nplot)/y0*abs(shr)) ! Have to take the complex conjugate to exactly reproduce the previous ! plots, which had an error because it didn't realize that ! the sin part = - imaginary part: yimag(l)=-yimag(l) enddo rf(ld)=r(ld-1) call c06fce(yreal,yimag,ld-1,work,ifail) ! fft ! shift origin to center of box do l=2,ld-1,2 yreal(l)=-yreal(l) yimag(l)=-yimag(l) enddo ! ! rotate phase at r=r0 phase=cmplx(yreal(1),yimag(1)) if(cabs(phase) > 1.e-30)then phase=phase/cabs(phase) else phase=(1.0,0.0) endif do l=1,ld-1 ylocal=cmplx(yreal(l),yimag(l))/phase yreal(l)=real(ylocal) yimag(l)=aimag(ylocal) enddo do l=1,ld-1 if (l <= ld/2) then ucft(l)=yreal(l+ld/2) usft(l)=yimag(l+ld/2) else ucft(l)=yreal(l-ld/2) usft(l)=yimag(l-ld/2) endif enddo ucft(ld)=0.0 usft(ld)=0.0 ! call hmplot2(ucft,usft,'Fourier transform of potential', ! & tim,ipage,mplot,nplot) ymin=0.0 ymax=0.0 call mnmx(ldb,ucft,ymin,ymax) call mnmx(ldb,usft,ymin,ymax) call setzer call plchlq(.5,.95,'Fourier transform of potential',20.,0.,0.) call plchlq(.5,.08,'k_par q R',15.,0.,0.) call mapsm(rf(1),rf(ldb),ymin,ymax) call curve(rf,ucft,ldb) call dashdb(21845) call curved(rf,usft,ldb) call finish(ipage, date) endif end subroutine hmplot2 subroutine hmplotc(a, title, tim, ipage) ! ! Plot mode amplitudes vs. ballooning angle for various m's. ! where a_ri is a set of complex mode coefficients. ! complex, dimension(:,:,:) :: a real, dimension(size(a, 1)) :: aloc real, dimension(2, size(a, 1), size(a, 2), size(a, 3)) :: a_ri real :: tim, ymin, ymax, xmin, xmax, yloc integer l, m, n, isympa, ipage character*70 title character*32 strtime character ctitle*80 integer ic, ltitle a_ri(1,:,:,:) = real(a(:,:,:)) a_ri(2,:,:,:) = aimag(a(:,:,:)) ! Loop over real and imaginary parts of the field: do ic=1,2 ymin=0. ymax=0. n=nplot2 do m=1,md do l=1,ldb yloc=a_ri(ic,l,m,n) if (yloc > ymax) ymax=yloc if (yloc < ymin) ymin=yloc enddo enddo xmin=r(1) xmax=r(ld) call setzer call plchlq(.5,.08,'ballooning coordinate',15.,0.,0.) if(ic == 1) then ctitle=title//' (real)' else ctitle=title//' (imag)' endif ltitle=index(ctitle,')') call plchlq(.5,.97,ctitle(1:ltitle),20.,0.,0.) write(strtime,'(a,i3,a,f7.2)') 'n = ',nplot2,', all m''s, time = ',tim call plchlq(.98,.93,strtime,15.,0.,1.) call mapsm(xmin,xmax,ymin,ymax) n=nplot2 do m=1,md do l=1,ldb aloc(l)=a_ri(ic,l,m,n) enddo call dashdb(dash(m)) call curved(r,aloc,ldb) enddo call finish(ipage, date) enddo ! end of do ic=2,1-1 end subroutine hmplotc subroutine hmplotcs(a, title, tim, ipage) ! ! Plot mode amplitudes vs. ballooning angle for various m's, ! where a_ri is a set of complex mode coefficients. For backwards ! compatibility, plot as "cosine" and "sine" parts, rather than ! as "real" and "imaginary" parts (including the sign difference). ! ! Eventually can drop this routine and use hmplotc instead. ! complex, dimension(:,:,:) :: a real, dimension(2, size(a, 1), size(a, 2), size(a, 3)) :: a_ri real, dimension(size(a, 1)) :: aloc real yloc real tim,ymin,ymax,xmin,xmax integer l,m,n,isympa,ipage character*40 title character*32 strtime character ctitle*80 integer ic,ltitle a_ri(1,:,:,:) = real(a(:,:,:)) a_ri(2,:,:,:) = aimag(a(:,:,:)) ! ! Loop over imaginary and real parts of the field: do ic=2,1,-1 ymin=0. ymax=0. n=nplot2 do m=1,md do l=1,ldb yloc=a_ri(ic,l,m,n) ! reverse imaginary part sign to reproduce previous plots exactly: if(ic == 2) yloc=-yloc if (yloc > ymax) ymax=yloc if (yloc < ymin) ymin=yloc enddo enddo xmin=r(1) xmax=r(ld) call setzer call plchlq(.5,.08,'ballooning coordinate',15.,0.,0.) ! put exactly the same title on as old cosine/sine version: if(ic == 2) then ctitle=title//' (sine)' else ctitle=title//' (cosine)' endif ltitle=index(ctitle,')') call plchlq(.5,.97,ctitle(1:ltitle),20.,0.,0.) write(strtime,'(a,i3,a,f7.2)') 'n = ',nplot2,', all m''s, time = ',tim call plchlq(.98,.93,strtime,15.,0.,1.) call mapsm(xmin,xmax,ymin,ymax) n=nplot2 do m=1,md do l=1,ldb aloc(l)=a_ri(ic,l,m,n) ! reverse imaginary part sign to reproduce previous plots exactly: if(ic == 2) aloc(l)=-aloc(l) enddo call dashdb(dash(m)) call curved(r,aloc,ldb) enddo call finish(ipage, date) enddo ! end of do ic=2,1-1 end subroutine hmplotcs subroutine hmplotk(ac, title, tim, ipage) ! ! Plot bounce-averaged complex field ac, vs. pitch angle kappa ! complex, dimension(:,:,:) :: ac character*70 title real tim integer ipage ! Used to be local variables, now in itg.cmn (be careful): ! real kap(kz) ! integer kd,kdpass ! local variables: real, dimension(:), allocatable :: ydata real yval, ymin, ymax, xmin, xmax, zkapmax integer ic, m, n, isympa, l character*25 strtime character ctitle*80 integer ltitle ! ! Define the kappa grid for the x-axis: ! kd=ldb/2 kdpass=ldb do l=1,kd kap(l)=sin(r(ldb/2+l)/2.) enddo if(epse > 0.0) then zkapmax=1./sqrt(2.*epse) else zkapmax=2.0 endif do l=kd+1,kdpass kap(l)=1.0+(1.0-kap(kd))+(l-kd-1)*(zkapmax-1.0)/(kdpass-kd) enddo ! define a mirror point on the other side of the maximum kappa ! to make differencing formulas easier: kap(kdpass+1)=zkapmax+(zkapmax-kap(kdpass)) allocate (ydata(kdpass)) ! Loop over real and imaginary parts: do ic=1,2 ! Find ymin and ymax for plot: ymin=0. ymax=0. ! do 10 n=1,nd n=nplot2 do m=1,md do l=1,kdpass if(ic == 1) yval=real(ac(l,m,n)) if(ic == 2) yval=aimag(ac(l,m,n)) if (yval > ymax) ymax=yval if (yval < ymin) ymin=yval enddo enddo ! ymax=1.1*ymax ! ymin=1.1*ymin xmin=kap(1) xmax=kap(kdpass) if (ymin == ymax) ymax=ymin+1. ! put exactly the same title on as old cosine/sine version: if(ic == 1) then ctitle=title//' (real)' else ctitle=title//' (imag)' endif ltitle=index(ctitle,')') call setzer call plchlq(.5,.08,'kappa (pitch angle)',15.,0.,0.) call plchlq(.5,.97,ctitle(1:ltitle),20.,0.,0.) write(strtime,'(a,i3,a,f7.2)') 'n = ',nplot2, ', time = ',tim call plchlq(.98,.93,strtime,15.,0.,1.) call mapsm(xmin,xmax,ymin,ymax) isympa=(ld-1)/10 ! plot only a single n (or plots for nonlinear runs take forever): ! do 12 n=1,nd n=nplot2 do m=1,md if(ic == 1) then do l=1,kdpass ydata(l)=real(ac(l,m,n)) enddo else do l=1,kdpass ydata(l)=aimag(ac(l,m,n)) enddo endif call dashdb(dash(m)) call curved(kap,ydata,kdpass) enddo call finish(ipage, date) enddo ! end ic=1,2 loop deallocate (ydata) end subroutine hmplotk subroutine cnplot0(acx,title,tim,ipage) ! ! Make contour plots of the field acx (complex mode coefficients). ! integer ntor,l1,l2 parameter(ntor=600) complex, dimension(:,:,:) :: acx real tim,ty(4*size(acx, 2)),tx(4*size(acx, 3)) real a2(4*size(acx, 3),4*size(acx, 2)),ator(ntor,ntor) real work(200*size(acx, 3)*size(acx, 2)), xx0, sc0, f1, f2, s1, s2 real rtor, xtor, ytor, thtor, q, alph, xx, zmax, zmin integer iwork(400*size(acx, 2)*size(acx, 3)) integer ipage, i, j, l, mr, nr, m, n, mpts, npts character*18 title character*30 maxstr logical annulus_data ! mpts=4*mz ! npts=4*nz mpts=4*md npts=4*nd do m=1,mpts ! ty(m)=2.0*pi*float(m-1)/float(mpts-1) ty(m)=2.0*pi*float(m-1-mpts/2)/float(mpts-1) enddo do n=1,npts ! tx(n)=2.0*pi*float(n-1)/float(npts-1) tx(n)=2.0*pi*float(n-1-npts/2)/float(npts-1) enddo ! ! contours at theta=0 ! l=(ld-1)/2 ! zmax=0.0 zmin=0.0 do i=1,npts do j=1,mpts a2(i,j)=0.0 do n=1,nd nr=nrr(n) do m=1,md mr=mrr(m,n) a2(i,j)=a2(i,j)+real(acx(l,m,n))*cos(-mr*ty(j)+nr*tx(i)) & -aimag(acx(l,m,n))*sin(-mr*ty(j)+nr*tx(i)) enddo enddo if (a2(i,j) > zmax) zmax=a2(i,j) if (a2(i,j) < zmin) zmin=a2(i,j) enddo enddo xx0=y0/(shr*z0)*float((nd-1)/2) sc0=xx0 if (y0 > sc0) sc0=y0 write(maxstr,101) zmin,zmax 101 format(' min: ',f7.2,' max: ',f7.2) ! Easy way ! call ezcntr(a2,npts,mpts) call cpseti('LLP',0) ! no labels on contours ! call cpseti('HIP',0) ! no high labels ! call cpseti('HLP',0) ! no low labels call cpseti('CLS - CONTOUR LEVEL SELECTION',0) call cpseti('NCL - NUMBER OF CONTOUR LEVELS',14) do j=1,14 call cpseti('PAI - PARAMETER ARRAY INDEX',j) xx=(zmax-zmin)/13.*float(j-1)+zmin call cpsetr('CLV - CONTOUR LEVEL VALUE',xx) call cpseti('CLU - CONTOUR LEVEL USE',3) if (xx < 0) call cpseti('CLD - CONTOUR LINE DASH PATTERN',3855) if (xx >= 0) call cpseti('CLD - CONTOUR LINE DASH PATTERN',65535) enddo call setzer call plchlq(.5,.95,title,20.,0.,0.) call plchlq(.65,.91,maxstr,9.,0.,-1.) ! call conop4('DAS=LSS','$$''''$$''''$''',0,0) ! call conop1('MES=OFF') ! call conop1('PMM=ON') ! call conop2('SSZ=ON',100) call cpseti('SET',0) call cpseti('MAP',3) call cpsetr('XC1',-xx0*pi) call cpsetr('XCM',xx0*pi) call cpsetr('YC1',-y0*pi) call cpsetr('YCN',y0*pi) call maps(-sc0*pi,sc0*pi,-sc0*pi,sc0*pi) call cprect(a2,npts,npts,mpts,work,200*mpts*npts,iwork,400*mpts*npts) call cpcldr(a2,work,iwork) ! call cplbdr(a2,work,iwork) call finish(ipage, date) if (title(1:9) == 'Potential') then write(9,*) '------ Contour data ' write(9,*) mpts,npts do j=1,mpts do i=1,npts write(9,*) a2(i,j) enddo enddo ! ! skip this section usually, it just makes fancy annulus data ! annulus_data=.false. if(.not. annulus_data) goto 33 do i=1,ntor do j=1,ntor ator(i,j)=0.0 xtor=float(i-1)/float(ntor-1)-0.500001 ytor=float(j-1)/float(ntor-1)-0.500001 rtor=sqrt(xtor**2+ytor**2)*2.*60. thtor=atan2(ytor,xtor) if ((rtor > 50.2).and.(rtor < 55.8)) then do l=1,ld-1 if ((r(l) < thtor).and.(r(l+1) > thtor)) then l1=l l2=l+1 endif enddo do n=1,nd nr=nrr(n) do m=1,md mr=mrr(m,n) q=qsf*(1.+(rtor-53.)/53.*shr) alph=-float(mr)*11.*q*(thtor-r0(m,n)) s1=thtor-r(l1) s2=r(l2)-thtor f1=s2/(s1+s2) f2=s1/(s1+s2) ator(i,j)=ator(i,j)+f1*( & real(acx(l1,m,n))*cos(alph)+aimag(acx(l1,m,n))*sin(alph) ) ator(i,j)=ator(i,j)+f2*( & real(acx(l2,m,n))*cos(alph)+aimag(acx(l2,m,n))*sin(alph) ) enddo enddo endif enddo enddo write(9,*) '------ more Contour data' write(9,*) ntor,ntor do i=1,ntor do j=1,ntor write(9,*) ator(i,j) enddo enddo 33 continue endif end subroutine cnplot0 subroutine cnplotk(a,title,tim,ipage,net,fave) ! ! k-space contour plot ! ! arguments: real, dimension(:,:,:) :: a character(*) :: title real tim,fave integer ipage,net ! local variables: integer l1,l2,n1 real, dimension(size(a, 2)) :: ty real, dimension(size(a, 3)) :: tx real, dimension(size(a, 3), size(a, 2)) :: a2 real, dimension(200*size(a, 2)*size(a, 3)) :: work real, dimension(:), allocatable :: tmpa integer, dimension(400*size(a, 2)*size(a, 3)) :: iwork real xx0,sc0,f1,f2,s1,s2 real ave,pfmax,pfmin,sdv,xx,zmax,zmin integer i,j,mr,nr,m,n,mpts,npts character*30 maxstr interface subroutine timeavp(seq,av,sdv,net,fave,smin,smax) real, dimension(:) :: seq real :: av, sdv, fave, smin, smax integer net end subroutine timeavp end interface mpts=md npts=nd ! mpts=mz ! npts=nz if (md <= 1) return if (nd <= 1) return ! ! time average ! zmin=0. zmax=0. allocate (tmpa(net)) do m=1,md do n=1,nd ! n1=nrr(n)+(nz-1)/2+1 n1=nrr(n)+(nd-1)/2+1 a2(n1,m)=0. do i=1,net tmpa(i)=a(i,m,n) enddo call timeavp(tmpa,ave,sdv,net,fave,pfmin,pfmax) a2(n1,m)=ave if (a2(n1,m) > zmax) zmax=a2(n1,m) if (a2(n1,m) < zmin) zmin=a2(n1,m) enddo enddo deallocate (tmpa) ! do m=md+1,mz ! do n=1,nz ! a2(n,m)=1.e-15 ! enddo ! enddo ! do m=1,mz ! do n=1,1+(nz-nd)/2 ! a2(n,m)=1.e-15 ! enddo ! enddo ! do m=1,mz ! do n=1+(nz+nd-2)/2,nz ! a2(n,m)=1.e-15 ! enddo ! enddo ! xx0=float(nz-1)*(abs(shr)*z0)/y0/float(nd-1) ! sc0=float(mz-1)/y0 xx0=abs(shr)*z0/y0 sc0=float(md-1)/y0 write(maxstr,101) zmin,zmax 101 format(' min: ',f7.2,' max: ',f7.2) call cpseti('LLP',0) ! no labels on contours call cpseti('CLS - CONTOUR LEVEL SELECTION',0) call cpseti('NCL - NUMBER OF CONTOUR LEVELS',14) zmax=zmax/3. do j=1,14 call cpseti('PAI - PARAMETER ARRAY INDEX',j) xx=(zmax-zmin)/13.*float(j-1)+zmin call cpsetr('CLV - CONTOUR LEVEL VALUE',xx) call cpseti('CLU - CONTOUR LEVEL USE',3) if (xx < 0) call cpseti('CLD - CONTOUR LINE DASH PATTERN',3855) if (xx >= 0) call cpseti('CLD - CONTOUR LINE DASH PATTERN',65535) enddo call setzer call plchlq(.5,.95,title,20.,0.,0.) call plchlq(.65,.91,maxstr,9.,0.,-1.) call cpseti('SET',0) call cpseti('MAP',3) call cpsetr('XC1',-xx0) call cpsetr('XCM',xx0) call cpsetr('YC1',0.) call cpsetr('YCN',sc0) call maps(-xx0,xx0,0.,sc0) call cprect(a2,npts,npts,mpts,work,200*mpts*npts,iwork,400*mpts*npts) call cpcldr(a2,work,iwork) ! call cplbdr(a2,work,iwork) call finish(ipage, date) if ((title(1:9) == 'Growth ra').or.(title(1:9) == 'Energy ').or. & (title(1:9) == 'Total dri')) then write(9,*) '------ k contour data ',title write(9,*) mpts,npts do j=1,mpts do i=1,npts write(9,*) a2(i,j) enddo enddo endif end subroutine cnplotk end module plots