pro script1, M, a_over_M, r, th, ph, t, dr, dth, dph, dl, nl, scene_format, pointset, cr, cg, cb, coo, coo_length, horizon, auto_radius, manual_radius, hcr, hcg, hcb, transp, make_r_stat, scr, scg, scb, out_file, append=append, verbose=verbose, labelx=labelx
    
    ; Bo Milvang-Jensen (milvang@astro.ku.dk)
    ; For the RelViz Web exhibition, June 1996
    ; Last rev.: 11-Jun-1996


    default, append,   0               ; Default is not to append
    default, verbose,  1               ; Default is to be verbose
    default, labelx,   0               ; Default is not to label the x axis

    ; The common block "my_common" is read by my version of pde
    common my_common, dummy_M, dummy_a

    COMMON hymanc, fold, first ; Used by hyman

    ; Calculate a. Note: 0 <= a_over_M <= 1, therefore 0 <= a <= M
    a = a_over_M * M

    ; Check if we are outside the horizon:
    r_H = M+(M^2-a^2)^0.5
    if (r lt r_H) then begin
      print, 'Error: you are inside the horizon; r = ', r, $
        ', r_H = ', r_H, ', returning.
      return
    end

    ; Check if # steps is too large (/tmp could otherwise be filled up)
    if (nl gt 10000) then begin
      print, 'Error: you have asked for too many steps (', nl, '), returning.'
      return
    end

    ; Put M and a in the common block "my_common"
    dummy_a = a
    dummy_M = M

    set_plot, 'null' ; No plots from IDL

    ; set dt to some value:
    dt = 1.0
    ; collect into fltarr(8)
    x  = [t, r, th, ph, dt, dr, dth, dph]

    ; We get a null geodesic by solving ds^2=0 for dt, which yields a
    ; 2nd order equation:
    ;
    kerr, x(1), x(2), x(3), g, cc, a=a, M=M ; calculate metric tensor g
    ;
    ; A plain 2nd order equation solution formula:
    aaa = g(0,0)
    bbb = 2*g(0,3)*dph
    ccc = g(1,1)*dr^2 + g(2,2)*dth^2 + g(3,3)*dph^2
    ;
    dt = (-bbb-sqrt(bbb^2-4*aaa*ccc))/(2*aaa)
    ;
    ds2 = g(0,0)*dt^2 + 2*g(0,3)*dt*dph + $
	  g(1,1)*dr^2 + g(2,2)*dth^2 + g(3,3)*dph^2

    if (verbose eq 1) then begin
      print, 'We calculate a dt of: ', dt
      print, 'This gives a ds^2 of: ', ds2, ' (should be 0)'
    end

    ; collect into fltarr(8)
    x  = [t, r, th, ph, dt, dr, dth, dph]

    ttt = fltarr(3,nl) ; ttt is used when calling writegeometry
    
    ; Make hyman work properly when script1 is called from call_script1
    ; (not a problem when called from the CGI script, since there one gets a
    ; new IDL all the time)
    first = 0.0

    ; step, and plot
    lambda=0.
    for il=0,nl-1 do begin
      hyman,x,lambda,dl
      ;
      ; See Schaum p. 50
      x_car = x(1)*sin(x(2))*cos(x(3)) ; x = r sin(theta) cos(phi)
      y_car = x(1)*sin(x(2))*sin(x(3)) ; y = r sin(theta) sin(phi)
      z_car = x(1)*cos(x(2))           ; z = r cos(theta)
      ;
      ; Save in variable we can use afterwards:
      ttt(0,il) = x_car
      ttt(1,il) = y_car
      ttt(2,il) = z_car
    end

    if (scene_format eq "Inventor") then vrml=0 else $
    if (scene_format eq "VRML"    ) then vrml=1 else begin
      print, 'Error: scene_format is neither "Inventor" nor "VRML", returning.'
      return
    end

    writegeometry, ttt, out_file, color=[cr, cg, cb], vrml=vrml, $
      pointset=pointset,append=append

    ; Add 3D widgets, that don't have their own script::
    openw, lun, out_file, /get_lun, /append

    ; Add coordinate system widget
    if (coo eq 'on') then begin
      printf, lun, 'Separator {'
      printf, lun, '        Scale { scaleFactor ', $
        coo_length, coo_length, coo_length, ' }'
      ;
      if (labelx ne 1) then begin
	coo_file = '/m/RelViz/WWW/milvang/Widgets_3D/coo_include'
      end else if (vrml eq 0) then begin
	coo_file = '/m/RelViz/WWW/milvang/Widgets_3D/coo_include_x_iv'
      end else begin
	coo_file = '/m/RelViz/WWW/milvang/Widgets_3D/coo_include_x_wrl'
      end	
      ;
      spawn, 'cat '+coo_file, coo_include
      printf, lun, coo_include, format='(a)'
      ;
      printf, lun, '}'
    end

    ; Add horizon widget:
    if (horizon eq 'on') then begin
      if (auto_radius eq 'on') then manual_radius = r_H
      printf, lun, 'Separator {'
      printf, lun, '        Material { diffuseColor ', $
        hcr, hcg, hcb, ' transparency ', transp, ' }'
      printf, lun, '        Sphere { radius ', manual_radius, ' }'
      printf, lun, '}'
    end

    ; Close file:
    close, lun
    free_lun, lun

    ; Add static limit / ergosphere widget:
    if (make_r_stat eq 'on') then begin
      write_r_stat, M, a, out_file, color=[scr,scg,scb], /append, nph=30, nth=15
    end


end