If I remeber my calculus the derivative of a series is not the same as the sum of the derivatives of the individual terms in the series i.e.

ds(y) = d(f1 + f2 + f3 + f4](y)
 dx                               dx

is not the same as

ds(y) = df1(y) + df1(y) + df1(y) + df1(y)
 dx       dx      dx        dx       dx

You need the first but what you are calculatibg is the second.

If you think about it the derivative (tangent) tn this case is the rate of change of y (height) at a position x (along width) and the easiest way of doing that is using a simple numerical approximation by looking at the difference in height between the point to the left and the point to the right.

I have modified your code to use this technique.
Note I have

• changed derivatives to tangents as it is more descriptive
• removed the horizon array because it is no longer needed
• used floats instead of doubles since the difference will not be noticeable and Processing methods all use floats.

1. final int w = 800;
2. final int h = 480;
3. float[] skyline;
4. PImage img;
5. int numOfDeriv = 800;
6. int derivModBy = 1; //Determines how many points will be checked
7. int time;
8. int timeDelay = 1000;
9. int iter;
10. float[] tangents;
11. 
    
12. public void setup() {
13.   noStroke();
14.   size(w, h);
15.   fill(0, 128, 255);
16.   rect(0, 0, w, h);
17.   terrain(w, h);
18.   fill(77, 0, 0);
19.   for (int i=0; i < w; i++) {
20.     rect(i, h, 1, -1*(int)skyline[i]);
21.   }
22.   time = millis();
23.   timeDelay = 100;
24.   iter =0;
25.   img = get();
26. }
27. 
    
28. public void draw() {
29.   if (iter == numOfDeriv) iter = 0;
30.   if (millis() > time + timeDelay) {
31.     image(img, 0, 0, width, height);
32.     strokeWeight(4);
33.     stroke(255, 0, 0);       
34.     point((float)iter*derivModBy, height-(float)skyline[iter*derivModBy]);
35.     strokeWeight(1);
36.     stroke(255, 255, 0);
37.     print("At x = ");
38.     print(iter);
39.     print(", y = ");
40.     print(skyline[iter]);
41.     print(", derivative = ");
42.     print((float)tangents[iter]);
43.     print('\n');
44.     lineAngle(iter, (int)(height-skyline[iter]), (float)tangents[iter], 100);
45.     lineAngle(iter, (int)(height-skyline[iter]), (float)tangents[iter], -100);
46.     stroke(126);
47.     time = millis();
48.     iter += 1;
49.   }
50. }
51. 
    
52. public void lineAngle(int x, int y, float angle, float length) {
53.   line(x, y, x+cos(angle)*length, y-sin(angle)*length);
54. }
55. 
    
56. public void terrain(int w, int h) {     
57.   //min and max bracket the freq's of the sin/cos series
58.   //The higher the max the hillier the environment
59.   int min = 1, max = 6;
60.   skyline =  new float[w];
61.   tangents = new float[w];
62.   //ratio of amplitude of screen height to landscape variation
63.   double r = (int) 2.0/5.0;
64.   //number of terms to be used in sine/cosine series
65.   int n = 4;
66.   int[] f = new int[n*2];
67.   //calculating omegas for sine series
68.   for (int i = 0; i < n*2 ; i ++) {
69.     f[i] = (int) random(max - min + 1) + min;
70.   }
71.   //amp is the amplitude of the series
72.   int amp =  (int) (r*h);
73.   for (int i = 0 ; i < w; i ++) {
74.     skyline[i] = 0;
75.     for (int j = 0; j < n; j++) {
76.       skyline[i] += ( sin( (f[j]*PI*i/h) ) +  cos(f[j+n]*PI*i/h) );
77.     }
78.     skyline[i] *= amp/(n*2);
79.     skyline[i] += (h/2);
80.   }
81.   for (int i = 1 ; i < w - 1; i ++) {
82.     tangents[i] = atan2(skyline[i+1] - skyline[i-1], 2);
83.   }
84.   tangents[0] = atan2(skyline[1] - skyline[0], 1);
85.   tangents[w-1] = atan2(skyline[w-2] - skyline[w-1], 1);
86. }
87. 
    
88. void reset() {
89.   time = millis();
90. }