module all_var 5
integer, parameter :: mrl=kind(1.0d0) !mrl=8
integer id,nextra,im,jm,icmax,iplot,iw,jw,ixy,iusexy,itrick
character*80 datafile,psfile,ctab
character*80 desc_file,ctab_file
double precision, allocatable :: x(:,:),y(:,:),q(:,:),extra(:),z(:,:),conlev(:)
integer, allocatable :: colia(:,:)
double precision header(6),magic_number
double precision xs,ys,smin,smax,qmin,qmax,qavg,qrms,&
xmin,xmax,ymin,ymax,xmn,xmx,ymn,ymx
character*6,allocatable :: cl(:)
character*5 adv
end module all_var
!***********************************************************
program d2ps ,25
use all_var
implicit double precision (a-h,o-z)
adv='yes'
magic_number=1.23d45
call read_input
call read_data
call make_eps_names
call open_ps_output
call diagnose_data
if (ixy.ge.0) call set_mx_mn
call set_scale
call set_size
call write_desc
if (iplot.eq.1.and.iusexy.ne.0)&
call pixel_header(int(512*xs+.9999),int(512*ys+.9999))
if (iplot.ge.2)&
call contour_header(int(512*xs+.9999),int(512*ys+.9999),iplot)
if (icmax.le.1.or.icmax.gt.1000) icmax=11
allocate (cl(0:icmax+3))
call colors(iwkid,ctab,iblk,icmax)
if (iplot.eq.1.and.iusexy.eq.0) then
if (im*jm.gt.10000) call big_im_plot
if (im*jm.le.10000) call im_plot
end if
if (iplot.eq.1.and.iusexy.eq.1) call pixel_plot
if (iplot.ge.2) then
nc=icmax-1 !number of contours
allocate(conlev(nc))
cint=(smax-smin)/nc ! the contour interval
ca=smin+cint/2. ! the first contour
do ii=1,nc
conlev(ii) = cint*real(ii-1)+ca
end do
write (*,*) ' xs,ys=',xs,ys
if (ixy.eq.0.and.iusexy.eq.1) then
call cplot(q,x,y,im,jm,conlev,nc,0,xs,ys,xmn,xmx,ymn,ymx,iplot)
else if (iusexy.eq.0.and.itrick.eq.1) then
call periodic_cplot(q,im,jm,conlev,nc,xs,ys,xmn,xmx,ymn,ymx,iplot)
else if (iusexy.eq.0.and.itrick.eq.0) then
call noxy_cplot(q,im,jm,conlev,nc,xs,ys,xmn,xmx,ymn,ymx,iplot)
else if (ixy.eq.1.and.iusexy.eq.1.and.itrick.eq.2) then
call b_cplot(q,im,jm,x,y,in,jn,conlev,nc,xs,ys,xmn,xmx,ymn,ymx,iplot)
else if (ixy.eq.1.and.iusexy.eq.1.and.itrick.eq.1) then
call periodic_xycplot(q,im,jm,x,y,in,jn,conlev,nc,xs,ys,xmn,xmx,ymn,ymx,iplot)
else
write (*,*) ' invalid options with contour plot? stopping'
stop
end if
end if
call close_ps_output
write (*,'(2a)',advance=adv) '..write ',trim(psfile)
deallocate(q,colia)
im=2*(icmax-1)
jm=1
allocate(q(im,1),colia(im,1))
do i=1,im
q(i,1)=smin+(smax-smin)*(i-.5)/(im)
end do
xs=1.
ys=.1
close (unit=97)
psfile=ctab_file
call open_ps_output
call im_plot
call close_ps_output
write (*,'(2a)') '..write ',trim(ctab_file)
contains
subroutine set_mx_mn 1
xmax=maxval(x)
xmin=minval(x)
ymax=maxval(y)
ymin=minval(y)
if (xmx-xmn.eq.0) then
xmx=xmax
xmn=xmin
end if
if (ymx-ymn.eq.0) then
ymx=ymax
ymn=ymin
end if
end subroutine set_mx_mn
subroutine read_input 1
open (unit=20,file='d2ps.input',status='old',err=100)
read (20,*) id,ixy,datafile,nextra,im,jm
read (20,*) xs,ys,xmn,xmx,ymn,ymx,ctab,icmax,smin,smax,psfile,iplot,iusexy,itrick
close (unit=20)
if (ixy.lt.0.and.iusexy.eq.1) then
write (*,*) ' ixy=',ixy,' iusexy=',iusexy,' stopping'
stop
end if
if (iplot.eq.1.and.iusexy.eq.1.and.ixy.eq.0) then
write (*,*) ' ixy=',ixy,' iusexy=',iusexy,' iplot=',iplot,' stopping'
stop
end if
return
100 write (*,*) ' no input file d2ps.input'
stop
end subroutine read_input
subroutine read_data 1
real headersp(6)
real, allocatable :: extrasp(:),qsp(:,:),xsp(:,:),ysp(:,:)
write (*,'(2a)',advance=adv)&
' d2ras will open ',trim(datafile)
if (id.eq.0) then
open (unit=19,file=trim(datafile),&
STATUS= 'old',form='formatted',err=100)
write (*,'(a)',advance=adv) '..opened formatted'
else
open (unit=19,file=trim(datafile),&
STATUS= 'old',form='unformatted',err=100)
write (*,'(a,i2)',advance=adv) '..opened unformatted, id=',id
end if
if (nextra.eq.0) then !look for header
if (id.eq.0) then
read (19,*) header(1)
else if (id.eq.1) then
read (19) headersp(1)
header(1)=headersp(1)
else if (id.eq.2) then
read (19) header(1)
else
write (*,*) ' id=',id,' should be 0<id<2'
stop
end if
write (*,*) ' header(1) was read'
if (abs(header(1)-magic_number).gt.1.d38) then
write (*,'(a,d15.5,d15.5,d15.5)',advance=adv)&
'..no header',magic_number,header(1),magic_number-header(1)
else
rewind (unit=19)
if (id.eq.0) then
read (19,*) header
else if (id.eq.1) then
write (*,'(a)',advance=adv) ' reading single precision'
read (19) headersp
header=headersp
else
read (19) header
end if
nextra=header(2)
im=header(3)
jm=header(4)
if (xs.le.0.0.or.xs.gt.1.0) xs=header(5)
if (ys.le.0.0.or.ys.gt.1.0) ys=header(6)
end if
end if
if (im.eq.0.or.jm.eq.0) then
write (*,*) 'headers says im=',im,' jm=',jm
stop
end if
write (*,'(a,3i4)',advance=adv)'..',nextra,im,jm
allocate(extra(nextra),q(im,jm),colia(im,jm))
if (id.eq.1) allocate(extrasp(nextra),qsp(im,jm))
if (ixy.ge.0) then
in=im+ixy
jn=jm+ixy
allocate(x(in,jn),y(in,jn))
if (id.eq.1) allocate(xsp(in,jn),ysp(in,jn))
end if
rewind(19)
if (id.eq.0) then
if (nextra.eq.0) read(19,*) q
if (nextra.ne.0) read(19,*) extra,q
if (ixy.ge.0) read (19,*) x,y
else if (id.eq.1) then
write (*,*) "will read single precision"
read(19) extrasp,qsp
extra=extrasp
q=qsp
deallocate(extrasp,qsp)
write (*,*) "did it ",extrasp
if (ixy.ge.0) then
read (19) xsp,ysp
x=xsp
y=ysp
deallocate(xsp,ysp)
end if
else
read(19) extra,q
if (ixy.ge.0) read (19) x,y
end if
close (unit=19)
write (*,*) "returning"
return
100 write (*,*) 'no file ',trim(datafile)
stop
end subroutine read_data
subroutine make_eps_names 1
idot=scan(psfile,'.',back=.true.)-1
if (idot.le.0) idot=len(trim(psfile))
desc_file=psfile(1:idot)//'_desc.eps'
ctab_file=psfile(1:idot)//'_ctab.eps'
psfile=psfile(1:idot)//'.eps'
end subroutine make_eps_names
subroutine open_ps_output 2
open (unit=97,file=trim(psfile),&
STATUS= 'unknown',form='formatted',err=100)
return
100 write (*,*) 'cannot open ',trim(psfile),' for output'
stop
end subroutine open_ps_output
subroutine close_ps_output 2
close (unit=97)
open (unit=94,file='catcom',status='unknown',form='formatted')
write (94,*) trim(psfile)
return
end subroutine close_ps_output
subroutine diagnose_data 1
qmin=minval(q)
qmax=maxval(q)
qavg=sum(q)/size(q)
qrms=sqrt(sum(q*q)/size(q)-qavg**2)
return
end subroutine diagnose_data
subroutine write_desc 1,18
character date_string*8,time_string*10
character*4 space
character*6,psc
call date_and_time(date_string,time_string)
space=achar(92)//'spc'
psc=') show'
open (unit=14,file=desc_file,status='unknown')
write (14,'(a)')&
'%!PS-Adobe-2.0 EPSF-2.0',&
'%%Creator: d2ps',&
'%%DocumentFonts: Helvetica',&
'%%BoundingBox: 0 0 512 96',&
'%%EndComments',&
'/Helvetica 12 selectfont'
iw=0
jw=0
call wmt
write (14,'(a,a,a)') '(datafile=',trim(datafile),psc
iw=iw+1;call wmt
write (14,'(a,a,a)') '(psfile=',trim(psfile),psc
iw=iw+1;call wmt
write (14,'(a,a,a)') '(ps_created=',&
trim(date_string)//' '//trim(time_string(1:4)),psc
iw=iw+1;call wmt
write (14,'(a,a,a)') '(ctab=',trim(ctab),psc
iw=iw+1;call wmt
write (14,'(a,es10.3,a)') '(xmin=',xmin,psc
call wmt
write (14,'(a,es10.3,a)') '(xmax=',xmax,psc
call wmt
write (14,'(a,es10.3,a)') '(ymin=',ymin,psc
call wmt
write (14,'(a,es10.3,a)') '(ymax=',ymax,psc
call wmt
write (14,'(a,es10.3,a)') '(xmn=',xmn,psc
call wmt
write (14,'(a,es10.3,a)') '(xmx=',xmx,psc
call wmt
write (14,'(a,es10.3,a)') '(ymn=',ymn,psc
call wmt
write (14,'(a,es10.3,a)') '(ymx=',ymx,psc
call wmt
write (14,'(a,es10.3,a)') '(qmin=',qmin,psc
call wmt
write (14,'(a,es10.3,a)') '(qmax=',qmax,psc
call wmt
write (14,'(a,es10.3,a)') '(qavg=',qavg,psc
call wmt
write (14,'(a,es10.3,a)') '(qrms=',qrms,psc
call wmt
write (14,'(a,es10.3,a)') '(smin=',smin,psc
iw=iw+2;call wmt
write (14,'(a,es10.3,a)') '(smax=',smax,psc
write (*,'(2a)',advance=adv) '..write ',trim(desc_file)
end subroutine write_desc
subroutine wmt 18
if (iw.gt.3) iw=0
if (iw.eq.0) jw=jw+1
write (14,'(i5,i5,a)')128*iw+2,98-16*jw,' moveto'
iw=iw+1
end subroutine wmt
subroutine set_scale 1
if (smin-smax.eq.0.) then
smin=qmin
smax=qmax
if (smin.lt.0.0.and.smax.gt.0.0) then
smax=max(abs(smin),smax)
smin=-smax
end if
end if
write (*,'(a,es9.2,es9.2)') '..',smin,smax
if (smin-smax.eq.0.) then
write (*,'(a)',advance=adv) '..no range'
stop
end if
end subroutine set_scale
subroutine set_size 1
if ((xs.eq.0.0.or.ys.eq.0.0).and.iusexy.eq.1) then
xs=maxval(x)-minval(x)
ys=maxval(y)-minval(y)
if (ys/xs.gt.0.3.and.ys/xs.lt.3.0) then
if (ys.gt.xs) then
xs=xs/ys
ys=1.
else
ys=ys/xs
xs=1.
end if
else
ys=1.
xs=1.
end if
end if
if (xs.le.0.0.or.xs.gt.1.0) xs=1.
if (ys.le.0.0.or.ys.gt.1.0) ys=1.
write (*,'(a,f6.3,f6.3)',advance=adv) '..scale ',xs,ys
end subroutine set_size
end program d2ps
subroutine big_im_plot 1,2
! too avoid memory overflows, writes out the image in stipes
use all_var
!call colors(iwkid,ctab,iblk,icmax)
call raster(q,colia,smin,smax,im,jm,icmax)
isize=512/im
jsize=512/jm
ib=im*isize
jb=jm*jsize
write(97,'(a)')&
'%!PS-Adobe-2.0 EPSF-2.0',&
'%%Creator: d2ps ',&
'%%DocumentFonts: Helvetica'
write (97,'(a)',advance='no') '%%BoundingBox: 0 0 '
write (97,*) ib,jb
write(97,'(a)')&
'%%EndComments',&
'%%EndProlog'
write (97,'(i5,1x,i5)',advance='no') isize,jsize
write (97,'(a)') ' scale'
do j=1,jm
write (97,'(i5,1x,i5)',advance='no') 0,min(1,j-1)
write (97,'(a)') ' translate'
write (97,'(a)') '<'
do i=1,im
write (97,'(a6)') cl(colia(i,j))
end do
write (97,'(a)') '>'
write (97,'(i5,1x,i5)',advance='no') im,1
write (97,'(a)')&
' 8',&
'[1 0 0 1 0 0]',&
'{}',&
'false',&
'3',&
'colorimage'
end do
write (97,'(i5,1x,i5)',advance='no') 0,-(jm-1)
write (97,'(a)') ' translate'
write (97,'(a)')&
'0 setlinewidth',&
'0 0'
write (97,*) im,jm
write (97,'(a)') '%rectstroke'
end
subroutine im_plot 2,2
use all_var
!call colors(iwkid,ctab,iblk,icmax)
call raster(q,colia,smin,smax,im,jm,icmax)
isize=xs*512/im
jsize=ys*512/jm
ib=im*isize
jb=jm*jsize
write(97,'(a)')&
'%!PS-Adobe-2.0 EPSF-2.0',&
'%%Creator: d2ps ',&
'%%DocumentFonts: Helvetica'
write (97,'(a)',advance='no') '%%BoundingBox: 0 0 '
write (97,*) ib,jb
write(97,'(a)')&
'%%EndComments',&
'%%EndProlog'
write (97,'(i5,1x,i5)',advance='no') isize,jsize
write (97,'(a)')&
' scale',&
'<'
do j=1,jm
do i=1,im
write (97,'(a6)') cl(colia(i,j))
end do
end do
write (97,'(a)')&
'>'
write (97,'(i5,1x,i5)',advance='no') im,jm
write (97,'(a)')&
' 8',&
'[1 0 0 1 0 0]',&
'{}',&
'false',&
'3',&
'colorimage'
write (97,'(a)')&
'0 setlinewidth',&
'0 0'
write (97,*) im,jm
write (97,'(a)') 'rectstroke'
end
subroutine pixel_plot 1,2
! (q,im,jm,smin,smax,xs,ys,ctab,icmax,&
! colia,psfile,x,y)
use all_var
!real*4 q(im,jm),x(im+1,jm+1),y(im+1,jm+1)
!integer colia(im,jm)
!character psfile*80,ctab*80,jnum*3
character jnum*3
!call colors(iwkid,ctab,iblk,icmax)
call raster(q,colia,smin,smax,im,jm,icmax)
xf=xs*12800/(xmx-xmn)
yf=ys*12800/(ymx-ymn)
do i=1,im
do j=1,jm
ibl=(x(i,j)-xmn)*xf
ibr=(x(i+1,j)-xmn)*xf
itl=(x(i,j+1)-xmn)*xf
itr=(x(i+1,j+1)-xmn)*xf
jbl=(y(i,j)-ymn)*yf
jbr=(y(i+1,j)-ymn)*yf
jtl=(y(i,j+1)-ymn)*yf
jtr=(y(i+1,j+1)-ymn)*yf
write(jnum,'(i3.3)') colia(i,j)
write (97,*) 'C'//jnum//' N'
write (97,*) ibl,jbl,' M'
write (97,*) ibr,jbr,' L'
write (97,*) itr,jtr,' L'
write (97,*) itl,jtl,' L P'
end do
end do
end
subroutine gscr(i,j,rv,gv,bv) 30,1
use all_var
character string*80,jnum*3,r*7,g*7,b*7
write(jnum,'(i3.3)') j
write(r(1:7),'(1x,f6.4)') rv
write(g(1:7),'(1x,f6.4)') gv
write(b(1:7),'(1x,f6.4)') bv
string='/C'//jnum//' {'//r//g//b//' C} bind def'
if (iusexy.gt.0.or.iplot.eq.2) write (97,*) trim(string)
ir=255*rv
ig=255*gv
ib=255*bv
write (cl(j),'(z2.2,z2.2,z2.2)') ir,ig,ib
return
end
subroutine hlsrgb(h,l,s,r,g,b) 3,3
real h,l,s,r,g,b,v1,v2
l=l/100.
s=s/100.
if (l.le.0.5) then
v2=l*(1+s)
else
v2=l+s-l*s
end if
v1=2*l-v2
if (s.eq.0.) then
r=l
g=l
b=l
else
r=F_valux(v1,v2,h+120.)
g=F_valux(v1,v2,h)
b=F_valux(v1,v2,h-120.)
end if
contains
function F_valux(v1,v2,hue) 3
real v1,v2,hue,F_valux
if (hue.gt.360.) hue=hue-360.
if (hue.lt.0.) hue=hue+360.
if (hue.lt.60.) then
F_valux=v1+(v2-v1)*hue/60.
else if (hue.lt.180.) then
F_valux=v2
else if (hue.lt.240.) then
F_valux=v1+(v2-v1)*(240.-hue)/60.
else
F_valux=v1
end if
return
end function F_valux
end subroutine hlsrgb
subroutine colors(iwkid,ctab,iblk,icmax) 1,33
implicit real (a-h,o-z)
character*80 ctab
if (iblk.eq.0) then
call gscr(iwkid,1,1.,1.,1.)
call gscr(iwkid,0,0.,0.,0.)
else
call gscr(iwkid,0,1.,1.,1.)
call gscr(iwkid,1,0.,0.,0.)
end if
call gscr(iwkid,2,.8,1.,1.)
call gscr(iwkid,icmax+3,1.,1.,.8)
if (trim(ctab).eq.'cyan_black_yellow') then
do k = 1,icmax
h = (k-1.)/(icmax-1.)
w=(1-h)*2
rv=max(0.,1.-w)
gv=.7*4*(h-.5)*(h-.5)+.3
bv=max(1.-(2*h),0.)
rv=max(rv,.7*bv)
call gscr(iwkid,k+2,rv,gv,bv)
end do
else if (trim(ctab).eq.'cyan_black_yellow_2') then
do k = 1,icmax
h = (k-1.)/(icmax-1.)
w=(1-h)*2
rv=max(0.,1.-w**2)
gv=.7*4*(h-.5)*(h-.5)+.3
bv=max(1.-(2*h)**2,0.)
rv=max(rv,.7*bv)
call gscr(iwkid,k+2,rv,gv,bv)
end do
else if (trim(ctab).eq.'cyan_black_yellow_3') then
do k = 1,icmax
h = (k-1.)/(icmax-1.)
w=(1-h)*2
rv=max(0.,1.-w**3)
gv=.7*4*(h-.5)*(h-.5)+.3
bv=max(1.-(2*h)**3,0.)
rv=max(rv,.7*bv)
call gscr(iwkid,k+2,rv,gv,bv)
end do
else if (trim(ctab).eq.'cyan_black_yellow_x') then
do k = 1,icmax
h = (k-1.)/(icmax-1.)
f=.9
w=(1-h)*2*f
g=2*h*f
rv=max(1.-w**3,0.)
gv=.7*4*(h-.5)*(h-.5)+.3
bv=max(1.-g**3,0.)
rv=max(rv,.7*bv)
call gscr(iwkid,k+2,rv,gv,bv)
end do
else if (trim(ctab).eq.'blue_grey_orange') then
do k = 1,icmax
h = (k-1.)/(icmax-1.)
rv=h
gv=.5
bv=1-h
call gscr(iwkid,k+2,rv,gv,bv)
end do
else if (trim(ctab).eq.'double_rainbow') then
do k = 1,icmax
h = real(k)/icmax
h=mod(2*h*360,360.)
call hlsrgb(h,50.,100.,rv,gv,bv)
call gscr(iwkid,k+2,rv,gv,bv)
end do
else if (trim(ctab).eq.'fruit4') then
do k = 1,icmax
h = (k-1.)/(icmax-1.)
rv=2*h-h*h
gv=.5*(1+cos(4*3.1415*h))*(1.-.25*sin(3.1415*h))
bv=1-h*h
call gscr(iwkid,k+2,rv,gv,bv)
end do
else if (trim(ctab).eq.'fruit3') then
do k = 1,icmax
h = (k-1.)/(icmax-1.)
g=h-.5
factor=.9999+.5*exp(-100*.25)-.5*exp(-100*g*g)
rv=(2*h-h*h)*factor
gv=.5*(1+cos(4*3.1415*h))*(1.-.25*sin(3.1415*h))*factor
bv=(1-h*h)*factor
call gscr(iwkid,k+2,rv,gv,bv)
end do
else if (trim(ctab).eq.'rainbow') then
do k = 1,icmax
h = real(k)/icmax
h=h*360+176.
call hlsrgb(h,50.,100.,rv,gv,bv)
call gscr(iwkid,k+2,rv,gv,bv)
end do
else if (trim(ctab).eq.'rainbow2') then
do k = 1,icmax
h = real(k)/icmax
g=0.5*(1-cos(3.141592*h))
r=g*360+176.
call hlsrgb(r,50.,100.,rv,gv,bv)
call gscr(iwkid,k+2,rv,gv,bv)
end do
else if (trim(ctab).eq.'fruit2') then
do k = 1,icmax
h = (k-1.)/(icmax-1.)
g=h-.5
factor=.9999+.5*exp(-200*.25)-.5*exp(-200*g*g)
rv=(2*h-h*h)*factor
gv=.5*(1+cos(4*3.1415*h))*(1.-.25*sin(3.1415*h))*factor
bv=(1-h*h)*factor
call gscr(iwkid,k+2,rv,gv,bv)
end do
else if (trim(ctab).eq.'fruit') then
do k = 1,icmax
h = (k-1.)/(icmax-1.)
g=h-.5
factor=.9999+.5*exp(-200*.25)-.5*exp(-200*g*g)
rv=(2*h-h*h)*factor
bv=.5*(1+cos(4*3.1415*h))*(1.-.25*sin(3.1415*h))*factor
gv=(1-h*h)*factor
call gscr(iwkid,k+2,rv,gv,bv)
end do
else if (trim(ctab).eq.'grey_tanh2') then
do k = 1,icmax
h = (k-1.)/(icmax-1.)
g=h-.5
tg=.25+.5*h+.2*(tanh(30*g)+tanh(15.))/(2*tanh(15.))
rv=tg
bv=tg
gv=tg
call gscr(iwkid,k+2,rv,gv,bv)
end do
else if (trim(ctab).eq.'grey_tanh') then
do k = 1,icmax
h = (k-1.)/(icmax-1.)
g=h-.5
tg=(tanh(8*g)+tanh(4.))/(2*tanh(4.))
rv=tg
bv=tg
gv=tg
call gscr(iwkid,k+2,rv,gv,bv)
end do
else if (trim(ctab).eq.'grey') then
do k = 1,icmax
h = (k-1.)/(icmax-1.)
rv=h
bv=h
gv=h
call gscr(iwkid,k+2,rv,gv,bv)
end do
else if (trim(ctab).eq.'grey_short') then
do k = 1,icmax
h = .3+.6*real(k-1)/(icmax-1.)
rv=h
bv=h
gv=h
call gscr(iwkid,k+2,rv,gv,bv)
end do
else if (trim(ctab).eq.'red_tanh') then
do k = 1,icmax
h = (k-1.)/(icmax-1.)
g=h-.5
tg=tanh(8*g)+tanh(4.)
factor=1./(2*tanh(4.))
!rv=tg*factor
rv=1.
bv=tg*factor
gv=tg*factor
call gscr(iwkid,k+2,rv,gv,bv)
end do
else if (trim(ctab).eq.'blue_white_red') then
do k = 1,icmax
h = (k-1.)/(icmax-1.)
green=min(1.9*h,1.9*(1-h))
red=.95
blue=.95
if (h.gt.0.5) blue=green
if (h.lt.0.5) red=green
call gscr(iwkid,k+2,red,green,blue)
end do
else if (trim(ctab).eq.'cyan_white_yellow') then
do k = 1,icmax
h = (k-1.)/(icmax-1.)
green=min(1.9*h,1.9*(1-h))
red=.95
blue=.95
if (h.gt.0.5) blue=green
if (h.lt.0.5) red=green
green=max(green,1-1.9*(1-h))
green=max(green,1-1.9*h)
call gscr(iwkid,k+2,red,green,blue)
end do
else if (trim(ctab).eq.'blue_black_red') then
do k = 1,icmax
h = (k-1.)/(icmax-1.)
green=0.
blue=max(1-1.3*h,0.)
red=max(1-1.3*(1.-h),0.)
call gscr(iwkid,k+2,red,green,blue)
end do
else if (trim(ctab).eq.'green') then
do k = 1,icmax
h = (k-1.)/(icmax-1.)
green=.18+.80*(-h*h+2*h)
red=.18+.80*h
blue=.18+.80*h*h
call gscr(iwkid,k+2,red,green,blue)
end do
else if (trim(ctab).eq.'red') then
do k = 1,icmax
h = real(k-1)/(icmax-1.)
red=.18+.80*(-h*h+2*h)
blue=.18+.80*h
green=.18+.80*h*h
call gscr(iwkid,k+2,red,green,blue)
end do
else if (trim(ctab).eq.'blue') then
do k = 1,icmax
h = real(k-1)/(icmax-1.)
blue=.18+.80*(-h*h+2*h)
green=.18+.80*h
red=.18+.80*h*h
call gscr(iwkid,k+2,red,green,blue)
end do
else !default
do k = 1,icmax
h = (k-1.)/(icmax-1.)
green=min(1.9*h,1.9*(1-h))
red=.95
blue=.95
if (h.gt.0.5) blue=green
if (h.lt.0.5) red=green
green=max(green,1-1.9*(1-h))
green=max(green,1-1.9*h)
call gscr(iwkid,k+2,red,green,blue)
end do
end if
return
end
subroutine raster(z,image,zmin,zmax,in,jn,icmax) 3
implicit double precision (a-h,o-z)
dimension z(in,jn)
integer image(in,jn)
zzmax=1.00000*zmax
fac=real(icmax-1)/(zzmax-zmin)
do j=1,jn
do i=1,in
ival=int((z(i,j)-zmin)*fac+3.5)
if (z(i,j).lt.zmin) ival=2
if (z(i,j).gt.zzmax) ival=icmax+3
image(i,j)=ival
end do
end do
return
end
subroutine pixel_header(ib,jb) 1
write (97,'(a)') &
'%!PS-Adobe-2.0 EPSF-2.0',&
'%%Creator: d2ps',&
'%%DocumentFonts: Helvetica'
write (97,'(a)',advance='no') '%%BoundingBox: 0 0 '
write (97,*) ib,jb
write (97,'(a)') &
'%%EndComments',&
'/M {moveto} bind def',&
'/L {lineto} bind def',&
'/R {rmoveto} bind def',&
'/V {rlineto} bind def',&
'/C {setrgbcolor} bind def',&
'/N {newpath} bind def',&
'/S {stroke} bind def',&
'% three choices follow: color and grid, color only, grid only',&
'/P {closepath gsave fill grestore 0 0 0 setrgbcolor stroke} bind def',&
'%/P {closepath gsave fill grestore} bind def',&
'%/P {closepath gsave grestore 0 0 0 setrgbcolor stroke} bind def',&
'%%EndProlog',&
'0 setlinewidth',&
'0 0 translate',&
'0.040 0.040 scale',&
'0 setgray'
return
end