This is a naïve transcription of the C code to OCaml. I found the graphics documentation for OCaml so simple that I stopped trying to make Java do the graphics. Here is a smattering of OCaml that might suffice. In place of the C array 
typedef struct{R x; R y; R xd; R yd;} pt;
pt a[m][n];
I adopted an array of arrays of records of four floats. Everything goes swimmingly until I get to
a.(i).(j).xd -:= p0 *. ddy;. I realized that OCaml lacks a feature that is often deprecated in C but is near the heart of the logic in my code. All arrays in OCaml are mutable but the fields of records are selectively mutable. I introduced the “+:=” operator because of the frequency of “+=” in the C code. The compiler then reminded me that “a.(i).(j).xd” evaluates to a float and not a reference. I can find neither a brief nor efficient workaround. The simplest I could find is: 
a.(i).(j+1).xd <- a.(i).(j+1).xd -. p1 *. dy.(j)
The text “a.(i).(j+1).xd” within “a.(i).(j+1).xd <-” has no meaning as an expression. The expression to the left of “+=” in C is an ‘l-value’ whose official definition I find obscure. The Algol 68 analog met this conundrum head on and the concepts were clear but not simple. Here is the bulkier slower working version whose output matches the C code exactly. (See the ap and av functions here. I have not decided to back port them.) The address of a record is calculated twice for each +:= and this dominates the cost of the calculation.

After adding graphics to this code (below) I begin to study the pragmatics of the version described just above.

The Movie

This for OCaml graphics on Mac OS X 10.7. This is a slightly interactive version. It opens an X11 window upon launch and waits for you to type in a decimal numeral followed by a CR. It will then compute and display that many time steps. If the value is 0 or there is no numeral that the previous count will be reused.

A still from the movie

This version maintains a positive pressure for triangles with negative area which is normally be avoided. Also added is an optional friction, “f”, on the boundary—a sort of viscosity.

This version is instrumented to learn about the distribution of pressures in order to color code triangles according to their pressure. Every time you type a “q” in place of a numeral, the pressures are ranked and 20 equal steps in the ranking are printed.

So far the programs are like the 650 version but with graphics added. Now we use new language features for a more modern version. This version shrinks by abstracting 2D vectors as tuples . We do arithmetic on the vectors with operators +^ -^ *^ *^. and !^. If x is a real and v a vector then “x *^ v” is scalar multiplication of the vector. “u *^. v” is the dot product between two vectors. The prefix operator “!^” rotates a vector counter-clockwise by 90 degrees. Vectors presumably halve the cost of computing addresses in the 2D array.

Here we reincorporate the graphics.

These modifications leave intact the extreme hacks to save space and work. We rely weakly on tail recursion only.

To be done: