PRO write_r_stat,M,a,file,color=color,append=append,vrml=vrml,indexedfaceset=indexedfaceset,nph=nph,nth=nth,transparency=transparency

; This procedure writes a geometry file showing the static limit
; It needs function r_stat

  ; Bo Milvang-Jensen (milvang@astro.ku.dk)
  ; For the RelViz Web exhibition, June 1996
  ; Last rev.: 07-Jun-1996


  default, color,    [1.0, 0.0, 1.0] ; Default color is purple
  default, append,   0               ; Default is not to append:
  default, vrml,     0               ; Default is to write Inventor, not VRML
  default, indexedfaceset, 0         ; Default is to write an "IndexedLineSet"
  default, nph,      30              ; number of points in phi direction
  default, nth,      15              ; number of points in theta direction
  default, transparency, 0.6         ; Default is 60% transparency

  openw,1,file,append=append ; Open file for writing, append if wanted:

  if (vrml ne 1) then begin
      header = [ $
      '#Inventor V1.0 ascii', $
      '' $
      ]
  endif else begin
      header = [ $
      '#VRML V1.0 ascii', $
      '' $
      ]
  endelse

  ; If we are not appending, print header:
  if (append ne 1) then printf,1,header,format='(a)'

  s1a = [ $
  'Separator {', $
  '    MaterialBinding {', $
  '	value	OVERALL', $
  '    }', $
  '    Material {', $
  '	ambientColor	0.2 0.2 0.2 ~' $
  ]

  printf,1,s1a,format='(a)'

  printf,1,'	diffuseColor	',string(color)
  printf,1,'	transparency	', transparency

  s1b = [ $
  '	specularColor	0.2 0.2 0.2 ~', $
  '	shininess	0 ~', $
  '    }', $
  '    Coordinate3 {', $
  '	point	[' $
  ]
  printf,1,s1b,format='(a)'

  f = r_stat(M,a,nph,nth)
  printf,1,f,form='(3f10.4,1h,)'
;  for i=0,nph-1 do begin
;    for j=0,nth-1 do begin
;      printf,1,f(*,j,i),form='(3f10.4,1h,)'
;    end
;  end

  s2= [ $
  '   ]' $
  ]
  printf,1,s2,format='(a)'

  s2= [ $
  '    }', $
  '    LightModel {', $
  '	model	BASE_COLOR', $
  '    }' $
  ]
  printf,1,s2,format='(a)'

  if (indexedfaceset ne 1) then begin
    s2a= [ $
    '    IndexedLineSet {' $
    ]
  end else begin
    s2a= [ $
    '    IndexedFaceSet {' $
    ]
  end
  printf,1,s2a,format='(a)'

  s2= [ $
  '        coordIndex [' $
  ]
  printf,1,s2,format='(a)'

  for i=0,nph-2 do begin   ; Note the -2 instead of -1 !
    for j=0,nth-2 do begin ; Note the -2 instead of -1 !
      printf, 1,  i   *nth+j  , ',', $
		  i   *nth+j+1, ',', $
		 (i+1)*nth+j+1, ',', $
		 (i+1)*nth+j  , ',',     ' -1,'
    end
  end

  s3= [ $
  '    ]', $
  '    }', $
  '}' $
  ]
  printf,1,s3,format='(a)'
  close,1

END