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