Previous Entry Next Entry→

July 20, 2005

Went to see the wonderful dentist Dr. Clapper yesterday.  The poor fellow thought he was just going to give me a crown and ended up doing a root canal to boot.  Yow!  Still, he is a paragon of doctorly expertise.  I am a much happier fellow for having a great dentist!

 The next exercise I’d like to tackle is to render a real mountain—say my big neighbor here, Mt. San Jacinto.  It looms large over every working day at College of the Desert.  I got the color relief map  of the Coachella Valley (see right) at the USGS web site, which also leads one to a custom data & download site where you can enter in your latitude and longitude intervals and get elevation data for that rectangular chunk.  I entered latitude bounds of 33.917 to 33.75 and longitude bounds from -117 to -116.5.  Now the web site says this about the data::   “Output data will be raw binary (with little-endian—Intel PC—byte order), is about a 180 by 360 degree area (21600 rows X 43200 columns = 933,120,000 grid cells), and each data value will occupy 2 bytes. The uncompressed data file size will be 1,866,240,000 bytes.”   So an interval of 33.917 to 33.75 degrees in latitude should pluck (33.917–33.75)/180*21600 ~ 20.04 rows of data (I got 21,  so that’s in the ball park) and (117-116.5)/360*43200 = 60, which is precisely what I got.  So each row is west to east data and the columns are configured north to south.  Here’s the Excel plot of these data:

To scale the clipping window, you need the maximum elevation for the vertical scale—which is easy enough to obtain, but for the horizontal scale, it might be useful to use the same unit of measure as the vertical scale.  Since we’re at an latitude of about 33.8 degrees, we are at a distance 6370*cos(33.8) = 5290km from the axis of rotation.  So the arclength along the circle of rotation is (117-116.5)*5290/360 = 7340 meters.

Well, this thing has eaten up a lot of time, so I don’t want to explain it in much detail other than to show the code and some cross-sections.

The code allows my to choose one of 21 different longitudinal (west-east) cross-sections such as the two below:

Here’s the code:

/* Draw a cross section of Mt. San Jacinto */

#include <GL/glut.h>   // This includes gl.h and glu.h

float west = 117.0, east = 116.5, north = 33.917, south = 33.75;

//DATA from usgs http://www.ngdc.noaa.gov/cgi-bin/seg/ff/nph-newform.pl/seg/topo/customdatacd

float mtSJ[21][60] = {{811,796,783,786,795,769,767,764,756,738,740,725,716,704,686,681,670,654,670,623,602,588,572,564,562,554,
542,522,506,490,473,456,438,428,417,406,398,380,373,362,351,354,368,374,377,364,344,315,297,291,259,242,228,228,230,226,225,224,222,219},

{792,786,766,797,783,799,759,763,758,740,728,715,700,690,673,680,670,655,675,670,714,617,562,545,536,522,510,499,484,475,461,447,440,511,445,
419,495,371,365,358,347,339,345,349,346,338,322,305,289,273,249,234,216,209,212,209,207,206,204,205},

{742,784,756,776,782,782,754,793,880,769,733,709,697,696,840,842,952,861,774,790,739,734,807,606,522,509,498,492,497,470,464,492,460,672,618,633,
794,502,361,350,344,345,332,327,324,322,301,288,277,256,252,229,213,201,193,212,212,193,189,190},

{715,698,759,789,757,753,738,730,743,813,776,772,772,829,971,1021,1050,980,921,948,953,975,886,843,730,663,717,938,747,555,542,760,777,833,1064,967,
672,454,382,394,381,389,540,464,523,296,290,279,267,232,243,231,211,205,193,174,179,180,175,176},

{676,664,717,728,687,747,729,681,750,732,742,780,898,1000,1119,1168,1183,1197,1135,1276,1215,1040,951,913,876,921,1058,1347,1047,923,718,960,1209,
1072,1233,935,867,753,492,463,649,629,606,771,549,335,301,284,281,248,225,227,207,193,184,168,160,163,163,162},

{669,655,698,758,691,657,701,661,693,711,803,892,945,1028,1219,1189,1190,1330,1292,1212,1190,1176,1049,1078,996,1032,1058,1097,1196,1180,1066,1019
,1309,1619,1294,1293,1208,996,616,532,817,1040,798,849,604,421,572,621,470,334,237,214,200,187,181,171,161,151,152,152},

{698,656,648,670,679,612,624,631,637,749,871,945,970,939,1015,1080,1156,1273,1306,1266,1226,1197,1206,1144,1090,1060,1086,1114,1248,1577,1457,1659,
1577,1932,1603,1382,1191,962,696,912,911,1280,1259,1179,715,586,717,830,653,510,337,217,201,186,176,166,158,149,142,161},

{682,585,597,603,607,601,638,649,644,675,819,792,875,943,1088,1165,1157,1168,1207,1214,1290,1258,1227,1214,1170,1162,1174,1216,1366,1574,1762,1902,
2061,1977,1649,1421,1278,899,1205,969,1065,1373,1500,1294,1200,1128,1041,1066,961,774,499,292,238,191,173,162,154,146,138,134},

{575,537,578,603,610,693,811,798,726,718,759,978,1065,1060,1090,1125,1083,1070,1074,1158,1266,1348,1307,1314,1237,1248,1312,1305,1398,1670,1843,2048,
2131,2032,1822,1734,1385,1437,1564,1236,1220,1570,1664,1624,1430,1350,1298,979,863,646,403,336,275,210,173,158,149,143,134,128},

{422,542,762,826,850,782,1013,993,856,805,848,998,1051,1071,1088,1139,1010,1093,1143,1206,1317,1508,1456,1441,1373,1346,1345,1475,1527,1804,2076,
2280,2210,2220,2128,1906,1558,1988,1963,1511,1633,1599,2040,2018,1634,1338,946,684,588,537,734,605,371,170,156,147,137,134,131,130},

{427,427,471,716,925,1089,1076,1038,977,987,1003,1026,1043,1071,1076,1124,989,1146,1202,1219,1324,1342,1446,1489,1521,1602,1598,1727,1933,2001,2178,
2245,2372,2493,2521,2111,1852,2253,2194,2097,2031,1717,2184,1997,1459,1314,906,1078,843,1027,820,566,248,155,139,133,130,126,124,121},

{436,433,434,477,675,804,894,1110,1039,1204,1085,1088,1077,964,1042,978,1045,1169,1153,1129,1169,1177,1289,1345,1479,1632,1754,1947,1968,2324,2174,
2006,2045,2021,2169,2486,2308,2599,2621,2567,2449,2056,2393,2058,1793,1427,1267,1421,1261,1280,968,732,525,371,129,124,119,117,114,115},

{430,433,447,443,440,478,668,911,991,100,928,845,96,825,883,874,958,1121,1191,1073,1076,1084,1222,1298,1319,1589,1721,1820,1876,2008,1965,1846,1840,
2008,2320,2444,2630,2999,3172,3119,2845,2704,2618,2493,2106,1946,1754,1497,1528,1361,1221,971,660,259,147,130,124,121,118,119},

{453,448,450,461,461,454,455,503,717,882,827,824,759,754,766,822,942,1129,1111,1053,983,1046,1109,1212,1331,1491,1743,1667,1804,1864,1793,1798,1859,
1939,2263,2364,2636,2869,3107,3006,2791,2702,2643,2569,2379,2106,1778,1362,1362,1170,1099,1035,585,228,153,129,123,117,113,105},

{456,452,457,465,469,465,463,468,469,564,615,613,628,646,635,659,788,933,1045,966,940,980,1153,1235,1328,1425,1617,1747,1611,1683,1567,1746,1749,2002,
2284,2605,2678,3000,3094,2945,2809,2729,2732,2670,2158,1786,1574,1417,1127,1034,1069,822,576,470,163,133,135,130,125,118},

{457,456,460,461,474,475,480,485,483,504,565,629,602,591,592,602,718,710,774,840,869,929,1036,1215,1227,1306,1434,1570,1402,1470,1677,1751,1979,2247,
2475,2508,2806,2909,2856,2633,2535,2462,2440,2334,2264,2045,1718,1408,1469,1336,1177,847,704,455,230,143,147,144,142,165},

{461,464,459,459,470,475,472,472,473,502,523,558,552,542,598,723,727,666,678,741,866,1045,1242,1337,1524,1620,1637,1359,1461,1749,1770,1834,1938,
2123,2166,2403,2463,2521,2404,2727,2563,2435,2296,2211,2138,2033,1880,1842,1619,1259,849,633,465,337,233,154,160,190,381,378},

{466,465,466,464,460,471,481,467,481,482,497,520,554,616,634,641,707,656,601,695,775,915,1268,1431,1529,1558,1488,1380,1446,1518,1580,1773,1896,
2040,2201,2294,2275,2175,2077,2462,2418,2397,2263,2287,2164,2096,2042,1632,1505,1439,1160,770,626,345,175,162,164,249,404,383},

{460,465,471,472,477,480,488,494,499,499,493,502,514,546,584,596,624,582,580,686,799,853,1090,1301,1393,1419,1264,1414,1496,1621,1670,1840,1970,
1950,1984,2032,1985,1889,2065,2532,2444,2426,2457,2377,1880,1741,1752,1691,1303,1450,1014,746,553,380,247,181,175,241,369,514},

{464,469,473,476,484,485,500,556,530,511,508,498,520,514,520,560,591,613,685,752,815,963,1177,1253,1248,1227,1180,1360,1578,1665,1782,1881,1806,
1793,1814,1866,1764,1941,2345,2517,2538,2513,2537,2329,2018,1660,1473,1318,1219,1185,1189,884,537,487,290,222,214,306,402,533}};

int maximum(float *a, int n) {

float max = a[0];

for (int i = 1; i <= n; ++i)

if (a[i] > max) max = a[i];

return max;

}

int longitude = 15;  // Change this for one of the other cross-sections: 0-20

int m = 60;

float maxElevation = maximum(mtSJ[longitude],m);

float longitudeLength = (west - east)*5290000.0/360.0;

float abscissa[60];

float increment = longitudeLength/60.;

void makeAbscissa() {

for (int i = 0; i<60; ++i)

abscissa[i] = i*increment;

}

void draw_axes() {

glColor3f(0.0,0.0,.0);

glBegin(GL_LINES);

glVertex2f(0.0,-100.);

glVertex2f(0.0,maxElevation);

glVertex2f(-100.,0.0);

glVertex2f(longitudeLength,0.);

glVertex2f(-50,maxElevation/2.);

glVertex2f(50,maxElevation/2.);

glVertex2f(longitudeLength/2.,-50.);

glVertex2f(longitudeLength/2.,50.);

glEnd();

glVertex2f(-40.,maxElevation);

glVertex2f(40.,maxElevation);

glVertex2f(0.,maxElevation+80.);

glVertex2f(longitudeLength,40.);

glVertex2f(longitudeLength,-40.);

glVertex2f(longitudeLength+80,0.);

glEnd();

}

void draw_segment(GLfloat ax, GLfloat ay, GLfloat bx, GLfloat by) {

glColor3f(0.0,0.0,1.0);

glBegin(GL_LINES);

glVertex2f(ax, ay);

glVertex2f(bx, by);

glEnd();

}

void draw_mountain() {

glColor3f(0.0,0.0,1.0);

for(int i = 0; i<60; ++i)

draw_segment(abscissa[i], mtSJ[longitude][i], abscissa[i+1], mtSJ[longitude][i+1]);

}

void display(void) {

glClear(GL_COLOR_BUFFER_BIT);

int n = 6;

glColor3f(1.0, 0.0, 0.0);

makeAbscissa();

glLineWidth(1.0);

draw_axes();

glRasterPos2f(100.,maxElevation-40.);// This labels vertical axis

glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_10, 'E');

glRasterPos2f(100.+120.,maxElevation-40.);

glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_10, 'l');

glRasterPos2f(100.+200.,maxElevation-40.);

glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_10, 'e');

glRasterPos2f(100.+300.,maxElevation-40.);

glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_10, 'v');

glRasterPos2f(100.+400.,maxElevation-40.);

glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_10, 'a');

glRasterPos2f(100.+500.,maxElevation-40.);

glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_10, 't');

glRasterPos2f(100.+600.,maxElevation-40.);

glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_10, 'i');

glRasterPos2f(100.+700.,maxElevation-40.);

glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_10, 'o');

glRasterPos2f(100.+800.,maxElevation-40.);

glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_10, 'n');

glLineWidth(3.0);

draw_mountain();

glFlush();

}

void init() {

//set clear color to white

glClearColor(1.0,1.0,1.0,1.0);

//set standard orthogonal (look straight at) clipping view

glMatrixMode(GL_PROJECTION);

// Scale window

gluOrtho2D(-100.0,longitudeLength+100.0,-100.0,maxElevation+100.);

}

int main(int argc, char** argv) {

//init mode and open window in upper-left corner

glutInit(&argc, argv);

glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);

glutInitWindowSize(400,400);

glutInitWindowPosition(0,0);

glutCreateWindow("dimple");

glutDisplayFunc(display);

init();

glutMainLoop();

}