program magmap c======================================================================= c UNIT 0: MAIN PROGRAM c======================================================================= c Driver program for making maps of magnetic fields or evaluating c field elements at specified points based upon spherical harmonic c models. Runs under the command language paradigm. c The program is built from 7 units: c UNIT 0: Main program c UNIT M: Mapping and plotting c UNIT I: Input commands c UNIT S: Spherical harmonic input and normalization c UNIT K: Decoding routines for Command Line Interface c UNIT Y: Generate and sum Ylm c UNIT P: Single point evaluation c parameter (inmx=200) parameter (mxdeg=96, mxb=(mxdeg+1)*(mxdeg+2)/2, mxa=mxb/10) complex blm,alm common /coeff/ lmax,lext,blm(mxb),alm(mxa) common /io/ inp,iout,iprint common /ndict/ iecho,nin,memory,istate(inmx) c ce 1000 call inpsh call getdat call point call mapper goto 1000 end c_______________________________________________________________________ subroutine projg(iproj, geo, u) c$$$$ calls nothing c Applies the forwd map projection from geographic coordinates geo c to the point u in R**2 on the map c dimension u(2), geo(2), x(3) common /orient/ center(2),O(3,3) data deg/0.01745329/,degby2/0.008726646/,pi/3.1415926/ ce if (geo(1) .gt. 90.0) then u(1)=999.0 return endif c cendg=center(2) c goto (1000, 2000, 3000, 4000, 5000), iproj c c Rectilinear projection: u(1)=longitude, u(2)=latitude c 0 <= u(1) <= 360; -90 <= u(2) <= +90 c 1000 u(2)=geo(1) u(1)=mod(geo(2) - cendg + 540.0, 360.0) - 180.0 return c c Lambert equal area: before rotations (0,0) in u plane is (0,0,1) in c R**3, and negative u(2) is the Greenwich meridian. The disk c u(1)**2 + u(2)**2 < 4 is mapped from the (punctured) sphere c Map onto unit sphere 2000 x(3)=sin(deg*geo(1)) rho =cos(deg*geo(1)) x(1)=rho*cos(deg*geo(2)) x(2)=rho*sin(deg*geo(2)) c Apply coordinate rotation O**T x if (O(3,3) .ge. 2.0) then x1=x(1) x2=x(2) x3=x(3) else x1 =O(1,1)*x(1) + O(2,1)*x(2) + O(3,1)*x(3) x2 =O(1,2)*x(1) + O(2,2)*x(2) + O(3,2)*x(3) x3 =O(1,3)*x(1) + O(2,3)*x(2) + O(3,3)*x(3) endif c c Map into Lambert plane: rho = 2*sin(theta/2) rho=sqrt(2.0*(1.0 -x3)/(x1**2 + x2**2)) u(1)= rho*x2 u(2)=-rho*x1 return c c Hammer-Aitoff equal-area projection: unit sphere is mapped on c the region u(1)**2/4 + u(2)**2 < 1 3000 clat=cos(deg*geo(1)) half=degby2*(mod(geo(2) - cendg + 540.0, 360.0) - 180.0) R=sqrt(1.0/(abs(1.0 + clat*cos(half)) + 1.0e-15)) u(1)=2.0*R*clat*sin(half) u(2)=R*sin(deg*geo(1)) return c c Orthogonal projection: before rotations (0,0) in u plane is (0,0,1) c in R**3, and negative u(2) is the Greenwich meridian. The disk c u(1)**2 + u(2)**2 < 1 is mapped into the (punctured) sphere. c Map onto unit sphere 4000 x(3)=sin(deg*geo(1)) rho =cos(deg*geo(1)) x(1)=rho*cos(deg*geo(2)) x(2)=rho*sin(deg*geo(2)) c Apply coordinate rotation O**T x if (O(3,3) .ge. 2.0) then x1=x(1) x2=x(2) x3=x(3) else x1 =O(1,1)*x(1) + O(2,1)*x(2) + O(3,1)*x(3) x2 =O(1,2)*x(1) + O(2,2)*x(2) + O(3,2)*x(3) x3 =O(1,3)*x(1) + O(2,3)*x(2) + O(3,3)*x(3) endif c Project upper hemisphere onto equatorial plane: rho = sin(theta) if (x3 .ge. 0.0) then u(1)= x2 u(2)=-x1 else u(1)=2.0 u(2)=2.0 endif return c c Mercator projection: u(1)=longitude (radians), u(2)=M c -pi <= u(1) <= pi; -pi/2 <= u(2) <= +pi/2 c 5000 u(2)=log(tan(degby2*(geo(1)+90.0))) u(1)=(mod(geo(2) - cendg + 540.0, 360.0) - 180.0)*deg if (abs(u(2)) .gt. pi) u(1)=999.0 return c end c_______________________________________________________________________ blockdata coast c Loads common /world/ with a coastline data set. Each real word c contains a lat-long pair recovered by: c lat=c(j)/10; long=1000*mod(abs(c(j),1.0) c Breaks signaled by lat=99 c common /world/ c(3217),nworld data nworld/3217/ c data (c(j),j=1,130)/94.28,90.2823,84.283,90.2837,95.2842, $100.284,110.285,115.287,119.288,110.288,121.29,115.291, $110.292,104.292,100.298,94.2989,86.2992,79.301,70.3016, $63.3026,59.3032,56.3064,51.3074,46.3081,42.3086,38.3092, $30.3093,20.3096,16.3102,8.31,0.3093,-9.3083,0.3091, $-5.3118,-10.314,-14.315,-20.3156,-27.316,-23.317,-29.318, $-33.3209,-40.3219,-46.3227,-51.3235,-59.325,-70.3252, $-80.3251,-89.3249,-99.3242,-107.323,-110.323,-114.323, $-120.322,-127.322,-132.321,-140.321,-150.321,-160.321, $-170.321,-177.321,-183.32,-191.321,-200.32,-210.319, $-220.319,-230.318,-234.315,-237.314,-245.313,-250.312, $-254.312,-263.311,-270.312,-280.311,-286.311,-292.31, $-300.31,-309.309,-317.309,-321.308,-328.308,-334.307, $-340.307,-347.306,-342.302,-349.302,-360.303,-370.303, $-376.303,-380.303,-384.302,-389.3,-394.298,-402.298, $-406.298,-411.297,-407.295,-412.295,-420.295,-427.296, $-436.295,-443.295,-449.295,-454.293,-460.292,-466.293, $-470.293,-476.294,-481.294,-484.293,-490.293,-496.292, $-501.292,990.37,-521.3,-516.299,-512.299,-516.302,-520.301, $-521.3,990.37,-501.292,-504.291,-510.291,-516.291,-522.292, $-530.292,-535.292,-540.293,-545.294/ data (c(j),j=131,260)/-550.294,-554.291,-551.29,-544.289, $-537.288,-530.287,-520.287,-512.286,-505.286,-497.286, $-489.286,-480.286,-470.286,-466.286,-460.286,-453.287, $-446.287,-440.287,-431.287,-421.288,-425.287,-433.286, $-427.286,-420.286,-413.286,-401.286,-391.287,-382.286, $-372.286,-360.287,-350.288,-340.288,-330.288,-320.288, $-310.288,-302.288,-296.289,-289.288,-283.289,-273.289, $-261.289,-250.289,-240.289,-230.289,-220.29,-210.29, $-200.29,-190.29,-180.289,-173.289,-165.288,-160.287, $-156.286,-152.285,-143.284,-135.284,-126.283,-117.283, $-109.282,-100.282,-90.2813,-80.2809,-70.2803,-64.2799, $-59.2788,-50.2788,-41.2787,-35.2795,-29.2801,-22.2791, $-10.2791,0.2798,10.2805,16.2812,27.2817,37.283,46.2828, $56.2827,66.2828,73.2822,83.2818,89.2813,84.28,990.37, $310.033,320.035,330.035,340.036,347.036,356.036,367.036, $361.034,360.033,990.37,350.033,345.033,350.034,355.035, $350.033,990.37,360.033,366.032,361.031,365.029,374.028, $380.027,388.027,394.026,402.027,410.03,420.034,405.038, $405.039,420.041,435.039,450.038,454.037,472.038,464.034, $452.036,450.034,457.033,461.033,466.031,464.03,435.029, $410.029,407.027,401.024,405.023/ data (c(j),j=261,390)/396.023,389.024,385.024,379.025, $372.024,373.023,990.37,450.047,464.049,479.05,479.052, $457.053,450.051,420.052,413.052,417.054,383.053,364.054, $361.052,376.049,398.05,435.047,450.047,990.37,356.024, $350.025,354.025,356.024,990.37,373.023,364.023,369.022, $366.022,370.022,379.021,383.021,390.021,399.02,406.02, $417.02,423.019,429.018,434.017,438.016,446.015,450.015, $456.014,450.013,443.012,439.013,435.014,429.014,423.015, $419.015,416.016,412.017,408.018,402.019,406.017,399.017, $395.017,390.017,383.017,379.016,387.016,990.37,381.016, $375.013,369.015,366.015,370.015,376.015,381.016,990.37, $387.016,397.016,405.015,410.014,414.013,420.012,427.011, $430.011,990.37,429.01,425.009,418.009,414.009,422.01, $429.01,990.37,412.009,407.01,400.01,391.01,400.009, $406.008,409.008,412.009,990.37,400.003,393.003,398.003, $400.003,990.37,430.011,440.01,444.009,438.008,431.007, $434.005,430.003,421.003,417.003,413.002,401,396.36,389, $383.36,377.359,369.358,364.355,361.355,990.37,104.107, $96.1061,90.1053,86.1049,94.1046,101.105,107.104/ data (c(j),j=391,520)/110.103,121.103,128.102,135.101, $120.1,110.099,100.099,92.0993,83.1001,74.1005,69.1009, $63.1024,53.103,43.1034,33.1033,23.1039,14.1041,20.1027, $25.1019,35.1012,43.1006,59.1004,70.0999,78.0993,82.0984, $91.0983,100.099,111.099,123.099,133.099,144.098,153.098, $163.098,170.097,163.096,990.37,135.093,125.093,117.093, $126.092,135.093,990.37,163.096,158.095,171.094,183.094, $190.094,196.093,202.093,211.092,220.092,230.091,220.09, $213.087,208.087,200.086,190.084,183.084,176.083,170.082, $164.082,159.081,153.08,144.08,134.08,126.08,117.08,106.08, $96.0789,90.078,82.0774,990.37,100.08,97.0805,89.081, $77.0816,69.0817,63.081,73.0796,87.0797,100.08,990.37, $82.0774,90.0764,102.076,113.076,122.075,133.075,144.074, $152.074,160.073,170.073,180.073,190.073,200.072,210.073, $217.072,224.072,214.072,210.071,220.069,225.069,230.07, $238.068,241.067,249.067,255.067,251.062,255.059,263.057, $270.057,264.055,271.053,277.052,286.051,294.051,301.05, $307.049,300.049,294.048,286.049,279.049,270.05,263.05, $256.05,253.051,261.051,256.052,250.052,254.052/ data (c(j),j=521,650)/241.052,250.055,257.056,263.057, $252.057,243.057,238.058,228.059,224.06,214.059,204.059, $190.058,180.057,170.055,163.053,156.053,152.052,149.051, $140.049,133.047,129.045,134.043,143.043,153.043,164.043, $173.042,180.042,192.041,201.04,208.039,217.039,225.039, $235.039,241.038,251.037,260.037,268.036,276.036,282.035, $288.033,300.033,990.37,300.033,290.033,280.034,270.034, $260.035,250.035,240.036,229.036,222.037,210.037,200.037, $189.038,183.038,174.039,163.039,156.04,150.041,146.041, $139.042,130.043,121.044,116.043,107.044,104.045,108.046, $111.048,116.05,120.051,109.051,100.051,90.0509,80.0502, $70.0497,60.0494,50.0486,40.0478,31.047,23.046,17.045, $10.044,0.0432,-10.0422,-20.0413,-26.0406,-37.04,-47.0395, $-55.0392,-62.039,-70.0398,-80.0397,-90.0398,-100.04, $-107.041,-118.041,-129.041,-140.041,-149.041,-156.041, $-163.04,990.37,-120.049,-128.05,-140.05,-150.05,-160.05, $-156.05,-165.05,-173.05,-184.049,-194.049,-206.049, $-215.048,-225.048,-235.048,-244.047,-252.047,-256.045, $-253.045,-248.044,-239.044,-229.043,-219.043,-212.044, $-206.044,-198.045,-187.044,-174.044,-163.045/ data (c(j),j=651,780)/ -156.047,-150.048,-136.048,-132.049, $-120.049,990.37,-163.04,-171.039,-179.037,-187.036, $-194.036,-200.035,-206.035,-214.035,-223.036,-230.036, $-240.036,-247.035,-254.033,-260.033,-270.033,-280.033, $-290.032,-300.031,-310.031,-316.03,-324.029,-332.028, $-337.027,-342.025,-349.02,-344.019,-338.018,-328.018, $-324.018,-315.018,-307.017,-297.017,-288.017,-281.016, $-272.015,-260.015,-250.015,-241.014,-230.014,-223.015, $-215.014,-206.013,-197.013,-188.012,-178.012,-167.012, $-154.012,-144.012,-135.013,-127.013,-124.014,-114.014, $-104.014,-98.0132,-90.0129,-82.0133,-70.0129,-61.0123, $-50.012,-40.0112,-30.0103,-20.0095,-8.0089,-1.0094,10.0098, $20.0101,31.0103,40.0097,46.009,50.0056,60.0051,64.004, $59.001,56,52.359,48.358,52.356,45.3529,51.3507,59.35, $66.349,71.348,80.3472,90.3469,100.346,110.345,120.345, $124.343,134.343,144.343,150.343,157.344,167.344,176.344, $188.344,196.344,203.344,211.343,221.343,230.344,240.344, $249.345,260.346,266.346,273.347,280.347,288.349,291.35, $296.35,302.351,308.35,316.35,321.351,329.351,333.351, $341.353,350.354,355.354,359.354,355.355,360.359,366.001, $369.003,373.01/ data (c(j),j=781,910)/ 366.011,362.011,356.011,350.011, $346.011,342.01,336.011,332.011,327.014,323.015,315.016, $310.018,306.019,303.019,310.02,320.02,327.021,323.023, $316.025,311.028,310.033,990.37,687.049,691.049,695.049, $692.05,688.05,687.049,990.37,800.048,804.047,809.051, $806.052,802.049,800.048,990.37,806.057,811.055,815.057, $811.059,807.058,806.057,990.37,803.062,808.06,803.062, $990.37,810.066,807.064,812.063,810.066,990.37,704.057, $708.054,712.054,721.052,726.053,731.053,737.054,742.055, $751.056,754.058,758.059,762.062,765.066,768.067,762.068, $758.066,755.063,751.062,745.059,740.058,734.057,730.056, $723.057,714.056,709.057,704.057,990.37,680.05,683.051, $689.054,683.054,680.053,675.053,678.054,682.055,685.057, $688.059,683.059,687.061,693.06,698.06,703.059,698.061, $693.064,690.066,686.067,681.068,687.069,690.068,693.068, $700.067,705.067,710.067,713.068,720.069,729.07,734.07, $730.07,726.073,717.072,709.073,699.073,689.073,684.074, $676.073,670.072,667.071,660.072,667.073,672.074,680.075, $685.074,688.075,700.074,706.074,713.073,719.073,725.075/ data (c(j),j=911,1040)/ 730.074,723.076,713.075,710.078, $713.077,719.076,723.077,720.081,717.081,723.082,730.081, $736.08,739.087,743.086,749.087,755.089,758.093,762.099, $764.1,990.37,788.098,794.094,798.094,802.092,805.092, $808.093,812.096,808.097,805.097,802.097,798.1,792.1, $788.099,788.098,990.37,779.099,784.1,789.101,793.102, $790.104,784.105,781.101,779.099,990.37,764.1,771.101, $774.102,771.104,767.107,763.113,759.114,753.114,749.113, $746.112,741.11,736.108,729.106,736.11,740.11,736.112, $741.113,745.112,741.113,736.113,732.119,728.123,732.123, $736.123,732.128,728.13,724.13,721.129,716.129,711.13, $707.131,715.132,719.133,715.133,719.14,724.14,990.37, $746.139,751.137,757.137,760.138,756.141,750.14,746.139, $990.37,750.144,753.142,758.143,767.142,759.143,756.145, $752.145,750.144,990.37,747.15,750.147,754.147,751.151, $747.15,990.37,732.144,736.141,739.142,734.144,732.144, $990.37,724.14,728.141,725.145,722.149,716.15,710.152, $707.159,703.16,698.159,694.161,699.168,694.169,691.168, $687.17,990.37,716.24,712.237,715.236,719.235,725.235/ data (c(j),j=1041,1170)/ 729.235,735.235,738.236,744.235, $741.243,737.244,733.243,729.242,726.241,719.239,716.24, $990.37,693.244,700.243,703.245,710.241,713.243,722.241, $728.242,732.244,726.245,730.247,727.249,723.248,727.25, $720.251,716.252,724.252,731.252,736.253,732.255,728.255, $724.255,719.255,715.256,710.255,706.256,701.259,696.259, $692.258,686.257,690.254,693.253,688.253,685.25,691.246, $693.244,990.37,688.243,683.246,679.245,676.247,679.25, $676.251,669.253,678.252,682.253,686.252,682.255,677.258, $683.261,687.261,693.262,698.262,695.263,690.264,686.265, $680.264,672.264,680.265,686.266,693.266,698.264,702.263, $706.264,713.264,718.265,724.265,733.264,740.265,732.269, $726.268,719.266,716.266,707.267,703.268,696.268,689.269, $682.27,990.37,717.264,712.262,718.26,722.259,725.258, $729.258,725.259,732.26,736.259,739.261,736.263,731.262, $726.263,717.264,990.37,746.265,750.263,756.265,752.267, $746.266,746.265,990.37,750.262,756.26,760.256,766.255, $758.258,764.258,757.262,750.262,990.37,750.254,746.249, $749.248,757.243,762.244,755.249,758.25,764.25,758.252/ data (c(j),j=1171,1300)/ 755.255,750.254,990.37,760.241, $757.24,761.237,764.238,770.24,773.241,770.244,765.244, $762.242,767.242,760.241,990.37,778.247,774.25,778.247, $990.37,782.25,785.249,782.25,990.37,785.255,789.256, $793.254,788.257,785.26,781.261,784.257,785.255,990.37, $780.262,788.262,785.265,782.265,778.264,780.262,990.37, $781.271,784.267,789.266,794.268,800.264,810.265,813.265, $808.269,805.27,797.273,793.275,788.273,781.271,990.37, $765.27,769.271,773.273,778.272,774.274,781.275,787.273, $794.275,798.273,803.274,807.281,814.27,818.269,822.273, $826.277,829.28,826.291,829.295,825.297,821.298,817.296, $813.295,809.293,805.291,800.289,796.288,792.284,780.284, $774.282,770.281,765.282,762.28,765.27,990.37,766.269, $762.271,756.271,750.28,745.28,754.268,763.267,767.263, $771.264,765.267,766.269,990.37,721.27,733.271,737.274, $732.275,735.276,728.284,725.284,720.286,716.286,708.29, $703.292,697.293,692.292,684.292,680.293,673.296,670.297, $665.299,660.298,654.297,647.296,653.295,659.294,664.293, $656.292,650.293,647.294,640.295,624.295,629.294,634.292/ data (c(j),j=1301,1430)/ 630.292,626.293,618.294,622.291, $627.29,636.288,640.287,645.286,641.283,648.282,653.283, $661.286,664.286,671.288,680.287,684.286,688.285,693.284, $702.282,698.281,702.273,707.271,710.272,721.27,990.37, $683.284,675.285,671.284,677.283,683.284,990.37,658.274, $654.275,650.277,647.277,639.279,636.28,640.277,635.276, $631.275,636.274,641.274,645.274,653.274,658.274,990.37, $-782.183,-778.2,-772.202,-768.209,-765.21,-762.211, $-756.216,-750.221,-745.224,-742.229,-738.236,-731.237, $-734.239,-738.247,-744.252,-738.253,-746.255,-750.259, $-743.259,-740.259,-735.259,-732.26,-728.258,-724.258, $-721.258,-717.26,-720.264,-726.265,-729.266,-726.268, $-731.271,-737.278,-731.28,-726.282,-729.285,-726.288, $-720.285,-717.284,-711.283,-714.285,-710.285,-706.286, $-700.285,-694.284,-691.287,-686.289,-689.29,-694.29, $-690.292,-683.293,-677.292,-670.291,-667.292,-670.293, $-666.294,-662.294,-654.296,-649.296,-646.298,-640.299, $-636.301,-630.302,-635.303,-639.301,-644.301,-648.301, $-654.301,-658.3,-663.3,-668.299,-674.299,-679.299,-684.299, $-689.299,-694.299,-699.299,-703.299,-712.3,-716.3,-721.3, $-725.3,-730.3,-734.3,-741.299/ data (c(j),j=1431,1560)/ -746.299,-751.299,-754.3,-761.303, $-764.305,-769.308,-773.311,-777.312,-780.315,-776.318, $-780.32,-777.324,-772.326,-766.33,-762.332,-757.333, $-751.335,-744.336,-740.338,-737.339,-730.34,-727.341, $-723.344,-720.345,-714.348,-710.349,-704.351,-700.358, $-690.36,-696,-699.001,-695.013,-691.015,-696.016,-699.02, $-696.029,-691.033,-686.034,-689.035,-695.037,-699.038, $-695.04,-689.04,-682.041,-677.044,-673.046,-670.049, $-663.05,-660.051,-657.053,-660.055,-665.057,-670.057, $-673.059,-677.066,-688.07,-684.072,-689.074,-693.077, $-689.078,-685.079,-681.079,-675.081,-671.082,-663.081, $-659.082,-664.083,-661.087,-665.089,-668.089,-665.092, $-661.095,-653.096,-650.096,-654.098,-658.099,-653.101, $-659.105,-664.108,-669.109,-664.11,-660.111,-657.112, $-660.115,-664.115,-668.117,-659.121,-663.123,-660.123, $-666.127,-670.128,-656.131,-659.131,-663.138,-667.14, $-671.143,-676.146,-679.147,-683.148,-690.156,-693.159, $-700.161,-704.161,-709.164,-706.167,-710.168,-717.17, $-723.171,-729.17,-734.169,-740.168,-743.167,-748.167, $-756.164,-766.164,-775.164,-783.165,-780.167,-776.167, $-779.178,990.37,682.27,691.271,687.272,676.272,671.273, $676.274,681.273,686.274,691.275/ data (c(j),j=1561,1690)/ 697.275,691.278,684.279,677.278, $667.278,662.276,665.275,662.274,654.273,650.273,641.272, $635.27,630.269,620.267,610.266,600.265,590.265,586.266, $577.267,569.268,565.272,559.273,553.275,540.278,529.278, $519.279,510.28,516.281,522.282,530.281,538.281,545.28, $553.283,558.283,566.284,574.283,580.282,584.281,591.282, $597.283,603.282,607.282,614.282,623.282,616.288,610.289, $600.29,591.291,587.291,581.292,584.293,589.294,598.295, $603.295,596.296,588.297,577.298,568.299,561.298,556.299, $551.3,547.302,539.302,535.303,524.304,518.304,513.303, $513.302,990.37,509.303,514.304,505.304,496.303,499.304, $495.304,484.306,476.306,480.307,472.307,465.307,476.306, $468.305,474.305,479.301,483.302,490.302,501.303,509.303, $990.37,498.296,494.298,490.299,498.296,990.37,513.302, $507.301,500.3,496.293,491.292,485.291,480.29,474.29, $480.291,485.292,489.293,479.295,469.295,460.296,456.297, $460.299,469.3,462.3,458.3,454.3,448.298,443.296,436.295, $443.294,449.295,456.296,450.294,444.292,437.291,428.289, $420.289,416.29,408.287,398.286,390.285,394.285/ data (c(j),j=1691,1820)/ 384.285,377.284,370.284,379.284, $383.284,374.284,370.284,366.284,358.284,351.284,346.283, $339.282,334.281,325.28,320.279,310.279,300.279,290.279, $284.28,274.28,266.28,258.28,250.28,255.279,266.278, $278.277,285.278,292.277,299.276,295.275,300.274,290.271, $297.268,294.268,290.265,286.264,282.263,270.263,260.263, $250.262,240.262,230.262,220.262,215.263,210.263,200.263, $190.264,182.266,187.268,196.27,209.27,212.271,206.273, $197.273,181.272,186.272,173.272,165.272,158.271,154.276, $150.277,140.277,130.277,120.277,111.276,100.277,92.278, $88.279,94.28,990.37,219.275,226.276,230.278,225.28, $217.282,212.284,206.285,200.286,205.283,212.282,219.28, $224.279,220.277,216.276,219.275,990.37,196.287,192.291, $183.292,174.288,180.288,184.286,192.287,196.287,990.37, $84.28,75.2803,78.279,82.2786,88.2764,92.2766,97.276, $100.275,110.274,116.274,122.273,127.273,132.272,139.27, $146.268,152.267,157.267,161.266,158.264,161.262,167.261, $173.259,176.259,181.258,187.256,193.256,197.255,205.254, $214.255,220.254,228.254,234.254,240.253,247.252,254.252/ data (c(j),j=1821,1950)/ 261.251,265.251,269.25,274.25, $281.249,286.248,293.248,300.247,310.247,314.247,317.246, $314.245,307.245,300.246,294.246,289.247,283.247,277.248, $271.248,266.248,260.249,250.249,245.249,236.25,231.25, $236.25,242.249,247.248,256.248,261.248,266.247,270.247, $280.245,287.246,293.245,296.245,305.244,313.244,321.243, $328.243,334.242,342.241,352.239,360.238,369.238,380.237, $385.237,390.236,398.236,406.236,414.236,422.236,428.235, $438.236,448.236,458.236,467.236,475.236,483.235,490.235, $494.234,499.233,507.232,513.232,519.232,527.231,534.23, $542.23,549.229,990.37,189.204,200.204,193.205,187.204, $189.204,990.37,540.228,531.227,525.228,520.229,528.228, $535.228,540.228,990.37,549.229,554.227,562.227,568.226, $576.225,582.224,586.222,592.221,597.22,602.216,607.214, $599.211,596.21,600.208,609.209,613.209,609.208,603.208, $597.207,590.206,589.207,990.37,584.207,580.207,575.205, $567.206,573.207,579.207,583.208,584.207,990.37,589.207, $582.206,578.205,574.204,570.203,566.202,559.201,556.2, $549.197,545.196,551.197,554.198,558.198,563.2,567.201/ data (c(j),j=1951,2080)/ 575.202,580.202,589.203,586.202, $589.202,584.201,590.2,585.198,593.198,602.198,596.197, $598.196,990.37,600.194,600.194,990.37,598.196,603.195, $609.195,614.194,618.194,624.195,629.195,633.196,638.199, $645.199,648.199,643.197,651.193,654.193,659.193,663.194, $660.198,663.199,667.197,671.197,675.196,678.195,683.193, $687.194,690.196,695.197,700.197,703.198,708.201,713.203, $708.204,712.205,708.205,705.208,702.213,698.218,694.221, $691.222,696.226,699.229,703.23,697.231,701.232,706.232, $701.233,697.233,701.235,696.236,692.241,688.243,990.37, $687.17,690.171,700.17,697.174,694.178,990.37,708.179, $713.179,708.179,990.37,694.178,689.18,990.37,686.181, $680.183,676.184,666.185,670.186,665.189,660.19,656.189, $649.188,646.187,641.187,644.186,648.184,653.184,662.181, $656.181,651.18,990.37,646.179,649.177,645.176,642.178, $631.179,626.18,622.176,617.174,610.172,605.171,599.17, $604.169,598.166,604.166,598.165,589.163,584.162,577.162, $573.163,566.163,559.163,563.163,554.162,547.162,540.16, $531.16,522.158,515.158,509.157,520.156,530.156,540.156/ data (c(j),j=2081,2210)/ 550.155,560.156,569.157,575.157, $579.157,584.159,589.16,595.16,599.161,604.162,608.164, $613.164,622.164,626.164,614.162,609.161,606.16,613.16, $619.16,616.159,619.159,614.157,607.157,603.155,595.154, $590.154,595.151,591.146,594.146,590.142,584.141,577.14, $571.139,563.138,557.137,551.136,547.135,538.137,542.138, $536.138,542.139,538.14,533.141,522.141,516.141,510.141, $500.141,490.14,990.37,490.142,500.142,511.142,521.142, $533.142,541.142,534.143,523.143,512.144,501.144,488.145, $493.144,484.143,476.143,470.143,463.144,467.143,460.142, $467.142,475.142,482.142,490.142,990.37,490.14,480.14, $470.139,460.138,452.137,443.136,435.135,428.134,433.132, $427.131,420.13,410.13,405.129,400.128,990.37,452.142, $448.143,443.143,436.145,429.145,421.143,426.142,418.141, $422.14,428.14,433.14,437.141,444.142,452.142,990.37, $415.141,407.142,400.142,390.142,384.142,381.141,373.141, $366.141,358.141,350.14,357.14,353.139,347.139,352.139, $347.138,343.137,339.136,334.136,340.135,347.135,339.132, $335.132,328.132,320.132,311.131,320.13,326.131,331.13/ data (c(j),j=2211,2340)/ 337.13,343.131,349.132,356.133, $360.136,369.137,374.137,367.137,371.138,378.139,385.14, $394.14,404.14,411.14,415.141,990.37,340.133,344.134, $338.135,333.134,327.133,334.132,340.133,990.37,400.128, $391.128,384.129,377.129,370.13,360.13,352.129,344.127, $353.126,359.127,368.126,378.126,388.125,394.125,398.125, $393.122,388.122,392.122,401.122,406.122,410.122,400.12, $392.119,386.118,380.118,374.119,370.12,378.121,374.122, $369.122,365.121,360.12,356.12,350.119,342.12,330.121, $319.122,309.122,304.121,300.122,290.122,281.121,270.12, $259.12,252.119,243.118,236.117,230.116,226.114,219.113, $214.111,205.111,210.11,990.37,200.11,209.111,193.111, $186.11,183.11,187.109,194.109,200.109,200.11,990.37, $210.11,217.109,210.107,200.106,190.106,180.106,173.107, $164.108,152.109,140.109,130.109,120.109,111.109,108.108, $104.107,990.37,220.121,229.12,237.12,243.12,251.121, $241.122,230.121,220.121,990.37,186.121,172.123,162.122, $154.121,141.122,129.124,134.123,139.123,131.123,139.122, $135.121,126.122,123.121,135.12,146.12,153.12,162.12/ data (c(j),j=2341,2470)/ 170.12,180.12,186.121,990.37, $125.124,113.126,104.125,100.125,96.124,92.123,98.1224, $106.123,119.122,115.123,110.123,119.124,125.124,990.37, $98.1255,89.1264,79.1267,70.1266,64.1263,71.126,62.126, $56.1255,60.1247,64.1241,70.1241,73.1244,79.1236,75.1235, $78.1227,69.122,80.1221,86.123,82.1242,86.1247,90.1249, $98.1255,990.37,-13.1307,-8.131,-4.1321,-8.134,-20.134, $-30.1349,-26.1363,-20.1371,-14.1378,-18.1387,-22.1398, $-26.1414,-30.1425,-34.1434,-38.1445,-43.1453,-48.146, $-54.1458,-60.1476,-64.148,-69.1472,-76.1477,-85.1484, $-90.149,-95.1496,-100.15,-106.15,-102.15,-94.147,-81.1461, $-77.1451,-83.1431,-90.1434,-94.1426,-90.141,-81.14, $-74.1383,-62.1385,-55.1382,-49.1373,-46.1356,-41.1347, $-34.1336,-40.1332,-34.1327,-28.1318,-24.133,-20.1337, $-16.1316,-13.1307,990.37,-55.1483,-50.1511,-41.1515, $-44.1523,-56.1519,-61.151,-55.1483,990.37,-28.1508, $-33.1522,-41.1532,-49.1529,-39.1524,-35.152,-28.1508, $990.37,-55.1548,-60.1554,-68.1558,-63.155,-55.1548,990.37, $-100.16,-94.1596,-100.161,-100.16,990.37,-203.164,-207.165, $-216.166,-223.167,-214.165,-203.164,990.37,-123.137, $-131.136,-143.136,-149.136/ data (c(j),j=2471,2600)/ -154.136,-160.137,-166.138,-170.139, $-174.14,-170.141,-160.141,-150.142,-140.142,-130.142, $-120.142,-110.142,-107.143,-117.143,-126.143,-133.144, $-143.144,-150.145,-160.145,-170.146,-180.146,-190.146, $-194.147,-200.148,-204.149,-211.149,-220.149,-224.151, $-233.151,-240.151,-247.152,-256.153,-264.153,-271.153, $-281.154,-290.154,-300.153,-310.153,-320.152,-328.152, $-333.151,-341.151,-350.151,-360.15,-369.15,-376.15, $-379.149,-386.147,-389.146,-384.145,-380.145,-384.144, $-388.144,-384.143,-380.14,-375.14,-370.14,-360.14,-354.139, $-348.139,-341.138,-351.138,-340.137,-334.138,-326.138, $-334.137,-340.136,-350.136,-337.135,-330.134,-323.134, $-315.131,-320.128,-326.125,-330.124,-337.124,-344.119, $-350.118,-341.115,-335.115,-331.116,-321.116,-314.116, $-303.115,-291.115,-280.114,-270.114,-260.113,-255.113, $-261.114,-250.114,-240.113,-234.114,-225.113,-219.114, $-223.114,-216.115,-210.116,-206.117,-202.118,-199.119, $-190.121,-180.122,-171.122,-163.123,-171.123,-160.123, $-153.125,-145.125,-139.126,-148.128,-138.13,-133.13, $-127.13,-124.131,-116.132,-112.132,-119.133,-125.136, $-120.136,-123.137,990.37,-408.145,-412.147,-419.148, $-428.148,-436.147,-429.146,-421.145,-415.145/ data (c(j),j=2601,2730)/ -408.145,990.37,-179.177,-172.178, $-179.179,-179.177,990.37,-345.173,-351.174,-359.175, $-369.175,-365.175,-376.176,-379.177,-374.178,-386.178, $-390.178,-397.177,-406.176,-411.176,-416.175,-407.175, $-401.175,-393.174,-380.175,-374.175,-367.174,-360.174, $-354.173,-345.173,990.37,-407.173,-412.173,-418.174, $-425.174,-432.173,-439.173,-450.171,-459.171,-465.169, $-461.167,-450.167,-440.168,-431.17,-426.171,-419.171, $-410.172,-407.173,990.37,0.128,10.1277,20.1282,10.1283, $6.1287,2.129,-9.1283,0.128,990.37,-32.128,-38.1308, $-32.128,990.37,-38.1266,-31.126,-36.1272,-38.1266,990.37, $-69.1341,-60.134,-55.1345,-62.1347,-69.1341,990.37, $-103.124,-93.124,-87.1252,-84.1271,-90.1261,-94.1254, $-101.125,-103.124,990.37,-97.119,-94.12,-101.121,-97.119, $990.37,-87.12,-83.1206,-80.123,-87.1228,-90.1217,-87.12, $990.37,-90.1159,-81.1161,-86.1191,-90.118,-90.1159,990.37, $-83.1156,-86.1143,-83.1131,-80.1106,-77.1096,-70.1063, $-67.1054,-60.106,-68.1085,-62.1107,-66.1117,-70.1124, $-77.113,-83.1156,990.37,0.1201,-8.1203,-13.1208,-8.1217, $-15.1223,-19.1218,-24.122,-31.1225,-37.1222,-41.123, $-46.1233,-53.1231,-46.1223,-49.122,-40.1216/ data (c(j),j=2731,2860)/ -35.121,-26.1211,-30.1204,-40.1204, $-47.1205,-55.1205,-46.1198,-35.1196,-28.1189,-20.1194, $-11.1195,-7.1199,0.12,10.1203,16.1251,5.1247,0.1201, $990.37,85.1173,94.1181,103.119,112.12,102.119,92.1185, $85.1173,990.37,0.1091,10.109,20.1095,16.1111,23.1113, $30.1118,40.1136,46.114,50.115,60.116,70.1169,65.1176, $60.1182,53.1192,50.1183,44.1186,35.1175,23.118,14.1184, $9.119,0.1176,-10.117,-16.1163,-21.1167,-31.116,-38.1154, $-34.1145,-30.113,-34.1117,-27.1116,-16.1099,-12.11,-4.1094, $0.1091,990.37,0.1035,-6.1032,-10.1042,-19.1047,-15.1054, $-30.1065,-40.1059,-50.1058,-58.1057,-49.1034,-40.1024, $-30.1015,-20.1009,-10.1004,0.0997,10.099,18.0989,23.0977, $31.0974,39.0967,44.0958,56.0951,50.098,40.0988,35.0995, $29.1002,21.1009,16.1019,11.1026,5.1036,0.1035,990.37, $361.355,364.354,370.354,380.351,385.351,394.351,402.351, $410.351,423.351,431.351,434.352,443.359,450.359,460.359, $469.358,474.358,485.355,496.358,493.36,990.37,513.35, $521.35,524.35,531.351,534.35,539.35,545.352,551.352, $546.355,542.355,532.354,526.354,521.354,517.352,513.35, $990.37,500.354/ data (c(j),j=2861,2990)/ 506.355,509.356,514.357,518.355, $528.356,532.356,540.357,544.356,549.357,546.355,550.355, $554.355,558.355,553.355,559.354,564.355,569.354,574.353, $571.354,577.354,585.355,579.356,576.356,571.358,566.358, $562.357,558.358,547.359,543.36,536,531,523.002,519.002, $516.001,512.002,508.001,504.358,500.355,500.354,990.37, $493.36,497,501.002,507.002,513.004,517.005,526.005, $521.005,526.006,533.006,538.009,547.009,554.009,565.008, $570.009,574.01,571.011,563.011,558.01,554.011,550.011, $557.011,549.012,545.012,551.01,544.01,539.011,546.017, $542.019,549.02,553.022,560.021,567.021,575.022,568.024, $572.025,582.024,586.023,591.024,594.025,599.029,604.029, $600.024,606.022,617.022,627.021,633.023,639.024,645.025, $651.025,656.025,651.022,646.021,637.021,629.019,623.018, $617.018,610.017,605.018,599.019,593.019,588.018,575.017, $990.37,579.019,573.019,569.019,576.019,579.019,990.37, $575.017,565.016,560.016,553.014,562.013,570.012,580.012, $591.011,595.011,589.01,583.009,579.008,587.006,596.006, $606.006,615.005,621.006,625.007,629.008,636.01/ data (c(j),j=2991,3120)/ 644.011,649.012,658.013,666.014, $671.015,678.016,683.016,689.016,686.017,689.018,694.018, $698.02,702.023,706.024,709.027,706.03,990.37,765.017, $770.015,774.014,780.014,783.016,787.012,791.011,796.011, $792.016,799.016,794.019,790.02,786.022,782.022,778.024, $774.024,778.022,784.02,780.019,774.018,765.017,990.37, $794.021,798.022,802.018,798.027,794.026,791.024,794.023, $794.021,990.37,706.03,700.03,697.031,693.033,689.036, $697.037,692.039,680.04,676.041,670.041,664.041,661.04, $664.035,657.035,648.035,644.035,639.037,644.037,648.037, $645.04,654.04,658.041,665.042,661.044,669.045,680.044, $686.043,678.047,674.045,670.046,667.046,676.048,680.05, $990.37,599.316,606.314,610.312,615.311,622.31,629.309, $636.309,641.308,651.308,656.307,659.306,664.307,670.306, $675.306,681.307,684.308,689.309,693.309,697.309,700.308, $697.308,693.307,701.305,707.305,710.308,713.308,717.307, $720.305,723.305,727.305,730.305,736.304,741.304,744.303, $748.303,754.301,758.301,762.298,766.291,769.289,773.291, $777.289,780.287,785.287,788.29,793.294,797.295,800.295/ data (c(j),j=3121,3217)/ 803.292,807.294,811.296,814.298, $819.299,816.301,819.301,823.305,820.308,825.308,822.312, $818.314,823.313,827.312,830.312,833.32,837.331,833.333, $830.336,827.339,823.338,820.328,817.326,820.329,817.333, $820.336,816.338,812.336,816.34,811.347,807.346,801.343, $796.342,791.341,781.341,775.341,772.342,768.342,762.338, $757.341,753.342,750.343,745.341,742.341,738.34,734.34, $729.338,724.338,721.338,724.336,720.337,716.338,704.338, $712.336,705.335,702.334,697.337,690.335,687.334,682.331, $672.327,667.326,662.325,655.323,649.32,642.32,636.32, $631.319,621.318,614.318,608.317,601.317,599.316,990.37, $633.341,637.339,646.338,650.338,653.338,658.336,663.337, $658.339,655.339,652.339,656.339,660.34,657.341,660.341, $665.344,662.345,657.346,651.346,646.346,641.345,637.343, $633.341,990.37/ end c_____________________________________________________________________ c======================================================================= c UNIT M: MAPPING AND PLOTTING c======================================================================= subroutine atlas(kind) c$$$$ calls getchr getarr getone rotate mult c Determines map projection (kind) and prepares various constants c for it in common /orient/. c Finds the array size and proportions and the appropriate step c size and limits for the projection; returned to caller in /scan/. c kind=1: Radian projection c kind=2: Lambert equal-area c kind=3: Hammer-Aitoff equal area (A) c kind=4: Orthographic projection c kind=5: Mercator projection c character*1 type dimension unit(3),axel(3),O1(3,3),O2(3,3),ulim(5) common /io/ inp,iout,iprint common /scan/ m,n,umin,vmin,du,factor,hash common /orient/ center(2),O(3,3) save type, em,iadd c data type /'R'/, em/40/, deg/0.01745329252/, axel(3)/0/ data unit /0,0,1/, enlarg/1/, pi/3.14159265/, iadd/0/ ce call getchr('projection', type, nf) kind=index('rlaomRLAOM', type) if (kind .eq. 0) then print '(3a)','>>> Unknown map projection: ',type, $ '. Magmap substitutes proj=r' type='r' kind=1 endif if (kind .gt. 5) kind=kind-5 c call getarr('center', center, 2, ncen) if (ncen .eq. 1) center(2)=center(1) c call getone('dimensions', em, ignor) c call getone('magnify', enlarg, ignor) factor=1.0/max(1.0, enlarg) if (enlarg .gt. 1.0 .and.iprint.ge.0) $print '(/a,f8.2)', $' Map enlarged by the factor ',enlarg hash=kind+1/(factor+100/(center(1)+100/(center(2)+.123))) c c Alternative specification: directly by a window call getarr('window', ulim, 5, jwin) if (jwin .eq. 5 .and. ulim(1).lt.ulim(2) .and. $ ulim(3).lt.ulim(4) .and. ulim(5).gt.0.0) then umin=ulim(1) vmin=ulim(3) du=ulim(5) m=1 + nint((ulim(4)-vmin)/du) n=1 + nint((ulim(2)-umin)/du) iadd=5 hash=-kind+1/(umin+1/(vmin+1/(du+.567))) elseif (jwin .eq. 0) then iadd=0 endif igo=kind + iadd c c Prepare the projections goto (1000, 2000, 3000, 4000, 5000, $ 1010, 2010, 3010, 4010, 5010), igo c c kind=1: radian projection 1000 m=nint(em) n=2*m - 1 umin=-180.0*factor vmin=-90.0*factor du=-umin/(m-1) c 1010 if (ncen .eq. 1) center(2)=center(1) if (iprint.ge.0) $print '(/a/a,f8.1)',' Rectilinear projection', $' Central meridian ',center(2) return c c kind=2: Lambert equal-area 2000 m=nint(em) n=m umin=-1.4142*factor vmin=umin du=-2.0*umin/(m-1) c 2010 if (center(1) .ge. 90.0) then O(3,3)=2.0 else c Build the rotation matrix used in projx; arrange meridian through c the map center to be vertical on the page; if center is one of the c poles, the Greenwich meridian is vertical axel(1)= sin(deg*center(2)) axel(2)=-cos(deg*center(2)) call rotate(center(2), unit, O1) call rotate(center(1)-90.0, axel, O2) call mult(3,3,3, O2,O1,O) endif c if (iprint.ge.0) $print '(/a/a,2f8.1)',' Lambert equal-area projection', $' Geographic coordinates of map center ',center c return c c kind=3: Hammer-Aitoff equal area 3000 m=nint(em) n=2*m - 1 umin=-2.0*factor vmin=-factor du=-umin/(m-1) O(1,1)= cos(deg*center(2)) O(2,2)=-sin(deg*center(2)) 3010 if (iprint.ge.0) $print '(/a/a,f8.1)',' Hammer-Aitoff equal-area projection', $' Central meridian ',center(2) return c c kind=4: Orthographic projection 4000 m=nint(em) n=m umin=-factor vmin=umin du=-2.0*umin/(m-1) c 4010 if (center(1) .ge. 90.0) then O(3,3)=2.0 else c Build the rotation matrix used in projx; arrange meridian through c the map center to be vertical on the page; if center is one of the c poles, the Greenwich meridian is vertical axel(1)= sin(deg*center(2)) axel(2)=-cos(deg*center(2)) call rotate(center(2), unit, O1) call rotate(center(1)-90.0, axel, O2) call mult(3,3,3, O2,O1,O) endif c if (iprint.ge.0) $print '(/a/a,2f8.1)',' Orthographic projection', $' Geographic coordinates of map center ',center return c c kind=5: Mercator projection 5000 m=nint(em) n=1.5*m - 1 du=2.0*pi*factor/(n-1) umin=-pi*factor vmin=-0.5*du*(m-1) c 5010 if (iprint.ge.0) $print '(/a/a,f8.1)',' Mercator projection', $' Central meridian ',center(2) return c end c_______________________________________________________________________ blockdata middle common /orient/ center(2),O(3,3) data center/90,0/, O(3,3)/2/ end c_______________________________________________________________________ subroutine projx(iproj, u, x) c$$$$ calls nothing c Applies the inverse map projection from a point u in R**2 on the map c plane to the unit vector in R**3 in the unit sphere dimension u(2), x(3) common /orient/ center(2),O(3,3) data deg/0.01745329/ c goto (1000, 2000, 3000, 4000, 5000), iproj c c Rectilinear projection: u(1)=longitude, u(2)=latitude c 0 <= u(1) <= 360; -90 <= u(2) <= +90 c 1000 cosla=cos(deg*u(2)) x(1)=cosla*cos(deg*(u(1)+center(2))) x(2)=cosla*sin(deg*(u(1)+center(2))) x(3)=sin(deg*u(2)) return c c Lambert equal area: before rotations (0,0) in u plane is (0,0,1) in c R**3, and negative u(2) is the Greenwich meridian. The disk c u(1)**2 + u(2)**2 < 4 is mapped into the (punctured) sphere 2000 sq=u(1)**2 + u(2)**2 R=sqrt(1.0 - 0.25*sq) x(1)=-R*u(2) x(2)= R*u(1) x(3)=1.0 - 0.5*sq if (O(3,3) .ge. 2.0) return c Apply coordinate rotation x1 =O(1,1)*x(1) + O(1,2)*x(2) + O(1,3)*x(3) x2 =O(2,1)*x(1) + O(2,2)*x(2) + O(2,3)*x(3) x(3)=O(3,1)*x(1) + O(3,2)*x(2) + O(3,3)*x(3) x(2)=x2 x(1)=x1 return c c Hammer-Aitoff equal-area projection: u(1)**2/4 + u(2)**2 < 1 is c mapped onto the unit sphere 3000 q=sqrt(2.0 - u(2)**2 - 0.25*u(1)**2) x(3)=q*u(2) cosla=sqrt(1.0 - x(3)**2) sinh=0.5*u(1)*q/(cosla + 1.0e-10) x2=2.0*sinh*sqrt(1.0 - min(1.0, sinh**2))*cosla x1=(1.0 - 2.0*sinh**2)*cosla c Apply coordinate rotation - only about z axis x(1)= x1*O(1,1) + x2*O(2,2) x(2)=-x1*O(2,2) + x2*O(1,1) return c c Orthogonal projection: before rotations (0,0) in u plane is (0,0,1) c in R**3, and negative u(2) is the Greenwich meridian. The disk c u(1)**2 + u(2)**2 < 1 is mapped into the (punctured) sphere 4000 R=sqrt(u(1)**2 + u(2)**2) if (R .gt. 1.0) then c (u, v) lies off the sphere; place image on limb x(1)=-u(2)/R x(2)= u(1)/R x(3)= 0.0 R=1.0 else x(1)=-u(2) x(2)= u(1) x(3)= sqrt(1.0 - R**2) endif if (O(3,3) .ge. 2.0) return c Apply coordinate rotation x1 =O(1,1)*x(1) + O(1,2)*x(2) + O(1,3)*x(3) x2 =O(2,1)*x(1) + O(2,2)*x(2) + O(2,3)*x(3) x(3)=O(3,1)*x(1) + O(3,2)*x(2) + O(3,3)*x(3) x(2)=x2 x(1)=x1 return c c Mercator projection: u(1)=longitude (radians), u(2)=M c -pi <= u(1) <= pi; -pi <= u(2) <= +pi c 5000 cosla=1.0/cosh(u(2)) x(1)=cosla*cos(u(1)+deg*center(2)) x(2)=cosla*sin(u(1)+deg*center(2)) x(3)=tanh(u(2)) return c end c_______________________________________________________________________ subroutine rotate(theta, en, a) c$$$$ calls nothing c Finds the 3 x 3 rotation matrix a that transforms vectors rotating c them about the axis en by an angle theta (degrees) anticlockwise c looking from origin along the axis. c c Rotation about the axis n, where n is unit by the angle A. c Written intelligibly the matrix is really c c R = I cos A + n n (1-cos A) + sin A n x c c where n x means "n cross". dimension en(3), a(3,3) data rad/0.01745329252/ c d=sqrt(en(1)**2 + en(2)**2 + en(3)**2) if (d .eq. 0) print *,'>>> Zero length rotation axis' c cost=cos (rad*theta) sint=sin (rad*theta)/d sinq=2.0*(sin(rad*theta/2.0)/d)**2 c a(1,1)=cost + sinq*en(1)*en(1) a(1,2)= sinq*en(1)*en(2) - sint*en(3) a(1,3)= sinq*en(1)*en(3) + sint*en(2) a(2,1)= sinq*en(2)*en(1) + sint*en(3) a(2,2)=cost + sinq*en(2)*en(2) a(2,3)= sinq*en(2)*en(3) - sint*en(1) a(3,1)= sinq*en(3)*en(1) - sint*en(2) a(3,2)= sinq*en(3)*en(2) + sint*en(1) a(3,3)=cost + sinq*en(3)*en(3) c return end c______________________________________________________________________ subroutine mult(m,n,k, a,b,c) c$$$$ calls nothing c Elementary matrix multiplication: c = a b where array sizes are c those appearing in the dimension statement below dimension a(m,n), b(n,k), c(m,k) c do 1300 i=1, m do 1200 j=1, k c(i,j)=0.0 do 1100 l=1, n c(i,j)=c(i,j) + a(i,l)*b(l,j) 1100 continue 1200 continue 1300 continue c return end c______________________________________________________________________ subroutine mapper c$$$$ calls atlas getchr getarr getone shsum projx plot c Calculates a field model on a regular array, which is in a specified c map plane. Writes the required field element to a diskfile parameter (pi=3.14159265358979, deg=180.0/pi) parameter (mxdeg=96, mxb=(mxdeg+1)*(mxdeg+2)/2, mxa=mxb/10) character*80 name, oldnam, code*1, type*1, nrm*1 complex blm,alm common /coeff/ lmax,lext,blm(mxb),alm(mxa) common /scan/ m,n,umin,vmin,du,factor,hash common /orient/ center(2),O(3,3) common /io/ inp,iout,iprint dimension u(2),r(3),fld(500) save code, shrnk, oldnam, kount, radius c data code/'R'/, bdum/1.0e8/, shrnk/1.0/, kount/0/ data type/'N'/, name,oldnam/'fort.9',' '/, radius/1.0/ c c Check to see if a map is needed ce call getchr('projection', type, nf) if (type .eq. 'N' .or. type.eq.'n') then print '(a)',' No map required in this run' return endif kount=kount + 1 c call getchr('normalize', nrm, ignor) call getone('igrf', yr, jgrf) if (nrm.eq.'S' .or. jgrf.gt. 0) then shrnk=1.0e-3 print '(/1x,a)', $ 'Gauss coeffs assmed to be in nT; then mapped fields are in muT' endif c c Get the file name for array call getchr('output', name, nnam) if (name .eq. oldnam) then print '(2a/a)', $'>>> Warning - Previous array of field values will be overwritten' $ ,'>>> A new filename should be set with an output command' elseif (nnam .le. 0) then if (kount .gt. 1) write(name,'(a5,i2)') 'fort.',88+kount print '(2a)','>>> No output file assigned; default to: ',name endif close (unit=9) open(unit=9, file=name) c c Get output array size and projection information; found in /scan/ call atlas(iproj) if (iprint.ge.0) $print '(/a,2i4)',' Dimension of output array',m,n c c Get the field element desired; map onto an index call getchr('field', code, nf) le=index('RXYZFHVID',code) if (le .eq. 0) then print '(3a)','>>> Unknown field element: ',code, $ '. Magmap assumes field=R' code='R' le=1 endif if (iprint.ge.0) $print '(/(a,a))',' Map of field element: ',code, $ ' Written to output file: ',name(1:35) if (iprint.ge.0) $print '(a,2f8.3,a,2f8.3,a,f8.3)',' Map plane = (', $umin,umin +(n-1)*du,')x(',vmin,vmin +(m-1)*du,'); dx=dy=',du c bmin=1e8 bmax=-bmin call getone('radius', radius, ignor) c c Scan through u plane: vmin,umin, etc come through /scan/ from atlas do 1500 i=1, m u(2)=vmin + (m-i)*du do 1400 j=1, n u(1)=umin + (j-1)*du call projx(iproj, u, r) ************************************************************ * call revers(iproj, u,r) ************************************************************ c Insert dummy value for illegal coordinates if (r(1)**2+r(2)**2+r(3)**2 .eq. 0.0) then fld(j)=bdum goto 1400 endif c Sum SH series r(1)=radius*r(1) r(2)=radius*r(2) r(3)=radius*r(3) c call shsum(r, lmax, blm, b1,b2,b3,V) if (lext.gt. 0) then call shsum(r, -lext, alm, exb1,exb2,exb3,exV) b1=b1 + exb1 b2=b2 + exb2 b3=b3 + exb3 V =V + exV endif c c Select the kind of field to be mapped goto (1010,1020,1030,1040,1050,1060,1070,1080,1090),le c Radial field desired: 'R' 1010 fld(j)=b1 goto 1400 c North field desired: 'X' 1020 fld(j)=-b2 goto 1400 c East field desired: 'Y' 1030 fld(j)=b3 goto 1400 c Downward field desired: 'Z' 1040 fld(j)=-b1 goto 1400 c Total field desired: 'T' 1050 fld(j)=sqrt(b1**2 + b2**2 + b3**2) goto 1400 c Horizontal field desired: 'H' 1060 fld(j)=sqrt(b2**2 + b3**2) goto 1400 c Scalar potential V/a: 'V' 1070 fld(j)=V goto 1400 c Inclination desired: 'I' 1080 fld(j)=atan(-b1/sqrt(b2**2+b3**2))*deg/shrnk goto 1400 c Declination desired: 'D' 1090 fld(j)=atan2(b3, -b2)*deg/shrnk c 1400 continue write(9,'(1p,g12.4)') (shrnk *fld(j),j=1,n) do 1450 j=1, n if (fld(j) .ne. bdum) then bmin=min(fld(j), bmin) bmax=max(fld(j), bmax) endif 1450 continue 1500 continue close(unit=9) c c Inform user about sizes if (iprint.ge.0) $print '(/3(a,1p,g12.3))', $' Field values between',bmin*shrnk,' and',bmax*shrnk,' muT' c c Display results if requested call plot(m,n, iproj, name) c return end c_______________________________________________________________________ subroutine launch(lk, script) c$$$$ calls system c CAUTION: the call to system is machine dependent and must be deleted c if taken out of the local SUN environment. c Runs the plotfile in script. character*64 script c if (script .eq. ' ') return close(unit=10) if (lk.eq.0) call system('sleep 2; contour <'//script) if (lk.eq.1) call system('sleep 2; color <'//script) c Any PostScript screen display can be substituted for pageview call system('pageview mypost &') return end c_______________________________________________________________________ subroutine getdat c$$$$ calls getchr getone getarr serial bang c Input routine. Accepts SH coeffs from a diskfile of the scalar c magnetic potential V. The default is fully normalized complex c values, integrating to unity over the sphere. The coeffs are stored c this way. These are the coeffs loaded into common /coeff/. c c Each input line should contain l, m, real, imag. c Or equivalently l, m, cosine coeff, sine coeff. c Negative m coeffs must not be provided: they are assumed to be c those appropriate to real omega. c The coeffs may be Schmidt normalized referred to the surface, c or one of 3 fully normalized conventions c c c parameter (mxdeg=96, mxb=(mxdeg+1)*(mxdeg+2)/2, mxa=mxb/10) character*80 name, style,line, norm*2 dimension sum(4),w(4) complex blm,alm common /coeff/ lmax,lext,blm(mxb),alm(mxa) common /io/ inp,iout,iprint common /earth/ rad(4) save name, norm, radius, year c data name/' '/, norm/' '/, radius/1/, pi/3.14159265/ data factor/1/, year/0/, minyr/1899/ c c c Check first if internally stored IGRF model is requested call getone('igrf', year, inunit) if (inunit.ge.1) then if (1900.0 .le.year .and. year.le. 2005.0) then inunit=nint(year) norm='S' if (iprint.ge.0) $ print '(a,i5)',' Field evaluated from IGRF-10 for ',inunit goto 1010 elseif (year .ne. 0.0) then print '(a,f7.1)','>>> Inappropriate year for IGRF model:', $ year stop endif endif c c Ascertain the name of the input file call getchr('file', name, nfile) if (nfile .le. 0 .and. name(1:1) .eq. ' ') then print '(a)',' ', $ '>>> Program cannot continue without filename for SH coeffs' stop elseif (name(1:1) .eq. '*') then if (iprint.ge.0) print '(/a)', $ ' SH coefficients input from the terminal' inunit=5 else if (iprint.ge.0) print '(/a,a)', $ ' SH coeffs from file: ',name(1:51) inunit=7 open (unit=7, status='OLD', err=2100, file=name) endif if (nfile .le. 0 .and. name(1:1) .ne. ' ') then if (iprint.ge.0) print '(a)', $ ' Program continues with SH coefficients previously input' return endif c call getchr('normalize', norm, ignor) if (norm .eq. ' ') then print '(a)', $ '>>> Program cannot continue without SH normalization' stop endif if (index('S F1F2G1U ',norm) .eq. 0) then print '(a)','>>> Normalization must be one of: S F1 F2 G1 U' stop endif if (iprint.ge.0) $print '(a,a)',' Normalization of coeffs: ',norm c 1010 call getone('radius', radius, nrad) if (radius .eq. 0) radius=0.54701 rad(4)=radius c c Set the earth dimensions: call getarr('earth', rad, 3, nrad) if (iprint.ge.0) print '(a,f9.3)', $' Equatorial radius (km) ',rad(1), $' Polar radius (km) ',rad(2), $' Ref radius aref (km) ',rad(3), $' Evaluation r/aref ',rad(4) c c Zero the arrays do 1100 n=1, mxb blm(n) = 0.0 if (n.le.mxa) alm(n)=0.0 1100 continue c c In contrast to other parameters, lmax is reset to maximum value c each time through the SH input process elmax=sqrt(2*mxb + 0.25) - 0.5 nexmx=sqrt(2*mxa + 0.25) - 0.5 call getone('lmax', elmax, nfound) lmax=elmax call getone('scale', factor, ignor) if (factor .ne. 1.0 .and.iprint.ge.0) $print '(a,g11.4)',' Coeffs scaled by',factor c c Check normalization f=sqrt(2.0*3.1415926535) rt2=sqrt(2.0) ff=factor c Norm is F2 if (norm .eq. 'F2') ff=factor*f*rt2 c c Check for serial data; read in file here call getchr('serial', style, nser) if (nser .ge. 0 .or. inunit .ge. minyr) then call serial(inunit, 0, ll, mm, glm, hlm) if (inunit .ge. minyr) style=' ' inunit=0 endif c call getone('list', dum, nlist) c lmx=0 kmx=0 bmx=0 if (nlist.ge.0 .and. inunit.ne.5 .and.iprint.ge.0) $print '(a)', ' ', ' Table of coefficients as read' if (inunit.eq.5 .and.iprint.ge.0) $print '(/(a))', ' Enter SH Coefficients NOW:', $' Enter l, m, glm, hlm each on a new line; l < -499 terminates' c c c Input the SH file one line at a time. do 1200 k=1, mxb c c Read the line of values keying on inuit for source if (inunit.eq. 7) then 1175 read (7,'(a80)', end=1250) line if (line(1:1) .eq. '%') then if (nlist.ge. 0 .and. iprint.ge. 0) print '(a)', line goto 1175 endif read (line, *, err=2000) ldeg,mm, glm, hlm elseif (inunit.eq. 5) then read (*,*, err=2000) ldeg,mm, glm, hlm elseif (inunit.eq. 0) then call serial(0, k, ll, mm, glm, hlm) ldeg=ll endif ll=abs(ldeg) c if (nlist.ge. 0 .and. iprint.ge. 0) $ print '(2i5,2f12.2)', ldeg,mm,glm,hlm if (ldeg.le. -500) goto 1250 if (ll.gt. lmax .or. -ldeg.gt.nexmx) goto 1200 lmx=max(lmx, ldeg) kmx=max(kmx, -ldeg) c Convert Schmidt semi-normalized to fully normalized coefficients if (norm .eq. 'S') then vl=f*(-1.0)**mm /sqrt(2.0*ll+1.0) glm = vl*glm hlm =-vl*hlm if (mm .eq. 0) glm =rt2*glm c Convert unnormalized to fully normalized coefficients elseif (norm .eq. 'U') then vl=f*(-1.0)**mm*sqrt(bang(ll+mm)/(bang(ll-mm)*(4.0*ll+2.0))) glm = vl*glm hlm =-vl*hlm if (mm .eq. 0) glm =2.0*glm c Norm is G1 elseif (norm .eq. 'G1') then vl=f*(-1.0)**mm glm = vl*glm hlm =-vl*hlm if (mm .eq. 0) glm =rt2*glm endif c Scale and store in temp array j=1 + mm + (ll*(ll+1))/2 if (ldeg.ge. 0) blm(j)=ff*cmplx(glm, hlm) if (ldeg.lt. 0) alm(j)=ff*cmplx(glm, hlm) bmx=max(abs(blm(j)), bmx) 1200 continue c c 1250 if (inunit.eq.7) close (unit=7) if (bmx .eq. 0) bmx=1 lmax=lmx lext=kmx if (iprint.ge.0) print '(a,i5)', $' Maximum degree of internal coefficients:',lmax, $' Max degree of external coeffs:',lext if (iprint.ge.0) $print '(a,f9.5)',' Normalized radius r/aref = ',radius c Check if dipole term(s) are to be dropped for this read only trm=3.0 call getone('nondipole', trm, nondip) if (nondip .ge. 0) then if (iprint.ge.0) then if (trm .eq. 3)print '(a)',' Dipole terms deleted if present' if (trm .ne. 3)print '(a)',' Axial dipole term deleted ' endif blm(2)=0.0 if (trm.eq. 3.0) blm(3)=0.0 endif c c List the fully normalized coefficients if requested call getone('list', dum, nlist) if (nlist .ge. 0) then if (iprint.ge.0) print '(/a,f9.5,a,f9.3,a)', $ ' Fully normalized SH coefficients referred to radius',radius, $ ' = ',radius*rad(3),' km' bscale=10.0**(-3*int((log10(bmx)-4.0)/3.0)) if (bscale .ne. 1.0 .and.iprint.ge.0) $ print '(a, 1p,g10.2,a)', $ '(Values in the table must be multiplied by ',1.0/bscale,')' do 1350 l=0, lmax j=1 + (l*(l+1))/2 if (iprint.ge.0) print 110, l,(blm(k)*bscale, k=j, j+l) 110 format(i4,6f12.2:/(f16.2,5f12.2:)) 1350 continue if (lext.gt. 0) then print '(a)',' External field coefficients' do 1360 l=1, lext j=1 + (l*(l+1))/2 if (iprint.ge.0) print 110, l,(alm(k)*bscale, k=j, j+l) 1360 continue endif endif c c Print weighted sums call getone('sums', dum, nesum) if (nesum .lt. 0) return j=0 do 1600 ns=1, 4 sum(ns)=0 1600 continue do 1800 l=1, lmax el=l w(1)=el+1 w(2)=(el+1)**2 w(3)=(el+1)*(2*el+1)**2*(2*el+3)/el w(4)=el*(el+1)**3 do 1700 m=0, l j=j + 1 bmagsq=2.0*abs(blm(j))**2 if (m .eq. 0) bmagsq=0.5*bmagsq do 1650 ns=1, 4 sum(ns)=sum(ns) + w(ns)*bmagsq 1650 continue 1700 continue 1800 continue c if (iprint.ge.0) print '(/a, 1p, g12.4, a)', $' Field energy:',sum(1), '*a**3/(2*mu0)', $' RMS radial field:',sqrt(sum(2)/(4*pi)), ' ', $' Lower bound on ohmic heating:',sum(3), '*a/(sigma*mu0**2)', $' RMS grad Br:', sqrt(sum(4)/(4*pi)),'/a' return c c Error return 2000 print '(a)','>>> Unable to read input line from disk' print '(a,2i5)','>>> Last successful l, m:',ll,mm stop 2100 print '(a)','>>> Cannot open the file ',name stop end c_______________________________________________________________________ function bang(k) c$$$$ calls nothing c Factorial of integer; no checks for negative argument or overflow bang=1.0 do 1100 j=1, k bang=bang * j 1100 continue return end c_______________________________________________________________________ blockdata tellus c Set default values in km for the equatorial, polar and reference c radii of the earth c WGS84 values below common /earth/ rad(4) data rad/ 6378.137, 6356.752, 6371.20, 1.0 / c end c_______________________________________________________________________ subroutine serial(inunit,k, l,m, br,bi) c$$$$ calls getchr getone c Reads serial data or returns values upon request. c If inunit = 0 returns kth pair (br, bi) in SH coefficient list and c If 0 < inunit < minyr reads from disk file to eof c If minyr < inunit uses inunit as year for IGRF field model c then corresponding degree and order are l, m. c c If nondipole field is specified, routine assumes series starts at l=2 c parameter (mxdeg=96, mxb=(mxdeg+1)*(mxdeg+2)/2, mxa=mxb/10) parameter (lmmx=2*mxb) character*20 style, line*80 common /igrf/ ibeg(23),ghlm(21*130+208) dimension b1(lmmx) save nondip, nozer, kmx, b1 data minyr/1899/ c c c Either fill b1 array or return the l,m pair if (inunit .gt. 0) then c c nondip=2 if dipoles are dropped, 0 otherwise dum=3.0 call getone('nondipole', dum, nondip) if (dum .ne. 3.0) print '(a)', $ '>>> All l=1 terms deleted in serial input mode' nondip=2 + max(-2, nondip) c c Either read data from a file or from igrf common if (inunit .lt. minyr) then c Reads coefficients serially, absent l, m labels, saved in b1 c until needed i1=1 do 1500 l=1, lmmx 1100 read (inunit,'(a)',end=1600) line if (line(1:1) .eq. '%') goto 1100 read (line, *, end=1200, err=5000) (b1(i),i=i1, lmmx) 1200 i1=i 1500 continue 1600 kmx=i1-1 c c Get igrf coefficients from common statement else k1=1 + (inunit - 1900)/5 ibgk=ibeg(k1) + 2*nondip kend= ibeg(k1+1)-1 do 2100 ik=ibgk, kend b1(ik - ibgk +1)=ghlm(ik) 2100 continue kmx=ibeg(k1+1) - ibgk endif c c Coefficients entered into b1 array if (inunit.eq.7) close (unit=inunit) call getchr('serial', style, ignor) nozer=min(1, index(style, 'out')) else c c c Retrieve SH pairs from b1. First find l, m; then work out index in b1 kk=k + nondip l=sqrt(0.251+2*kk) - 0.5 m=kk - (l*(l+1))/2 n=2*(kk - nondip) - 1 - nozer*(l - nondip/2) if (m .eq. 0) n=n+nozer if (n .gt. kmx) l=-999 br=b1(n) bi=b1(n+1)*min(1,m) endif return c c Error return 5000 print '(a,i5)','>>> Error in serial data on line ', l print '(a)',line stop end c_______________________________________________________________________ subroutine inpsh c$$$$ calls clear icheck common /io/ inp,iout,iprint c Input routine for commands to read SH coeffs and plot maps c c Reads from the standard input until eof or code "continue" or "quit". c Saves lines in the input store /store/ for later retrieval by c getarr, getone or getchr. c c Prints a glossary upon request c parameter (inmx=200) character *80 input(inmx),line,cmd*4 common /store/ input common /ndict/ iecho,nin,memory,istate(inmx) c write(*,'(a)') ' ', $' =================',' ', $'Enter commands for spherical harmonic coefficients (? for help)' c do 1500 l=nin+1, inmx read(*,'(80a)', end=3000) line if (line(1:5) .eq. 'conti' .or. line(1:4) .eq. 'exec') goto 2000 c List a glossary of codes if (line(1:1) .eq. '?') write(*, '(a/a/(2x,a))') $'Enter commands from the following list:', $' (M) means mandatory information',' ', $'?: Remind me of the command list again', $'execute: Quit reading commands and begin calculations', $'quit: End calculations', $'file filename: Input the filename of coeffs', $'serial [without zeros] SH coeffs in natural order without l,m', $'normalization X: X=S Schmidt, X=F1 fully normalized, integral=1', $' X=F2 integral=4*pi, X=G1 Gravity normal, U=unnormalized (M)', $'igrf year: Use internal IGRF coeffs in place of file input', $'lmax: Accept coefficients up to this degree', $'nondipole: Delete dipole terms from all subsequent input', $'list: List SH file as read and in F1 normalization', $'radius r: Calculate field & plot map on sphere radius = r*aref', $'earth a,b,aref: Equatorial, polar, ref radius (km)', $' ---------', $'output filename: Diskfile to receive map array', $'silent: Suppress terminal output henceforth (except error mssg)', $'dimensions m, n: Array dimensions' if (line(1:1) .eq. '?') write(*, '(2x,a)') $'center lat, long: Geographic coordinates of map center', $'projection p: p=l Lambert equal area, p=r radian, p=a Aitoff, ', $' p=o orthographic, m=Mercator ', $'magnification f: Draw map with scale f times bigger than normal', $'coast [file]: Lat-long pairs for coast; blank =internal dataset', $'symbols file type ht [color=X]: File of lat-longs for symbols', $'field E: E=one of (X, Y, Z, H, F, R, V/a, I, D)', $'scale s: Multiply SH coeffs by s before plotting', $'site lat, long, alt: Evaluate field elements at site with given', $' geographic coords and altitude in meters', $'sums: Print certain weighted sums of coeffs', $'contour [file]: Prepare a contour script in file and run it', $'color [file]: Prepare a color contour script in file and run it', $'clear cmd: Remove all instances of cmd from stack', $' ' if (line(1:1) .ne. ' ') then if (line(1:4) .eq. 'quit') goto 3000 if (line(1:4) .eq. 'echo') iecho=-iecho if (line(1:4) .eq. 'sile') iprint=-1 c Alternative file designators suppress each other if (line(1:4) .eq. 'file') call clear('igrf') if (line(1:4) .eq. 'igrf') call clear('file') c Alternative plotting formats suppress each other if (line(1:4) .eq. 'cont') call clear('colo') if (line(1:4) .eq. 'colo') call clear('cont') c nin=nin + 1 call icheck(line) input(nin)=line istate(nin)=0 endif 1500 continue write(*,'(a,i5,a/a)')'>>> Command stack full: only',inmx, $' allowed','>>> Recompile MAGMAP with a larger inmx' stop c 2000 if (iprint.ge.0) then write(*, '(5x,a)') ' ','===================', $ (input(i)(1:75), i=1, nin), '===================',' ' else write(*,'(/a/a)')' Magmap screen output suppressed', $ '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~',' ' endif call getchr('clear', cmd, klr) c Remove cleared commands if (klr .gt. 0) call clear(cmd) if (klr .gt. 0) call clear('clear') c return c 3000 print '(/a)',' Magmap run complete' call plot(0,0,0,' ') stop end c______________________________________________________________________ subroutine icheck(line) c$$$$ calls nothing c Checks the command string line (1:4) against a list; appends a c warning if command is not present in the list character*80 line, com*4 ce com=line(1:4) if (0 .eq. @+index('cent clea coas colo cont dime eart fiel file igrf', com) @+index('list lmax magn nond norm outp plot proj radi scal', com) @+index('seri site sums symb wind', com) $) line=line(1:20)//' <<<<<<< Unrecognized command' return end c______________________________________________________________________ blockdata iounit parameter (inmx=200) common /io/ inp,iout,iprint common /ndict/ iecho,nin,memory,istate(inmx) data inp/5/, iout/6/, iprint /1/, iecho/-1/, nin/0/ end c______________________________________________________________________ blockdata igref c IAGA IGRF-10: 1900-2005. All degree 10, except 2005=deg 13 c See common /igrf/ ibeg(23), $y1900(130), y1905(130), y1910(130), y1915(130), y1920(130), $y1925(130), y1930(130), y1935(130), y1940(130), y1945(130), $y1950(130), y1955(130), y1960(130), y1965(130), y1970(130), $y1975(130), y1980(130), y1985(130), y1990(130), y1995(130), $y2000(130), y2005(208) data ibeg/1,131,261,391,521,651,781,911,1041,1171,1301, $1431,1561,1691,1821,1951,2081,2211,2341,2471,2601,2731,2939/ c data y1900/ $-31543,0,-2298,5922,-677,0,2905,-1061,924,1121, $1022,0,-1469,-330,1256,3,572,523,876,0,628,195, $660,-69,-361,-210,134,-75,-184,0,328,-210,264,53, $5,-33,-86,-124,-16,3,63,0,61,-9,-11,83,-217,2, $-58,-35,59,36,-90,-69,70,0,-55,-45,0,-13,34, $-10,-41,-1,-21,28,18,-12,6,-22,11,0,8,8,-4, $-14,-9,7,1,-13,2,5,-9,16,5,-5,8,-18,8,0,10, $-20,1,14,-11,5,12,-3,1,-2,-2,8,2,10,-1,-2,-1, $2,-3,0,-4,2,2,1,-5,2,-2,6,6,-4,4,0,0,-2,2, $4,2,0,0,-6/ data y1905/ $-31464,0,-2298,5909,-728,0,2928,-1086,1041,1065, $1037,0,-1494,-357,1239,34,635,480,880,0,643,203, $653,-77,-380,-201,146,-65,-192,0,328,-193,259,56, $-1,-32,-93,-125,-26,11,62,0,60,-7,-11,86,-221, $4,-57,-32,57,32,-92,-67,70,0,-54,-46,0,-14,33, $-11,-41,0,-20,28,18,-12,6,-22,11,0,8,8,-4,-15, $-9,7,1,-13,2,5,-8,16,5,-5,8,-18,8,0,10,-20, $1,14,-11,5,12,-3,1,-2,-2,8,2,10,0,-2,-1,2, $-3,0,-4,2,2,1,-5,2,-2,6,6,-4,4,0,0,-2,2,4, $2,0,0,-6/ data y1910/ $-31354,0,-2297,5898,-769,0,2948,-1128,1176,1000, $1058,0,-1524,-389,1223,62,705,425,884,0,660,211, $644,-90,-400,-189,160,-55,-201,0,327,-172,253,57, $-9,-33,-102,-126,-38,21,62,0,58,-5,-11,89,-224, $5,-54,-29,54,28,-95,-65,71,0,-54,-47,1,-14,32, $-12,-40,1,-19,28,18,-13,6,-22,11,0,8,8,-4,-15, $-9,6,1,-13,2,5,-8,16,5,-5,8,-18,8,0,10,-20, $1,14,-11,5,12,-3,1,-2,-2,8,2,10,0,-2,-1,2, $-3,0,-4,2,2,1,-5,2,-2,6,6,-4,4,0,0,-2,2,4, $2,0,0,-6/ data y1915/ $-31212,0,-2306,5875,-802,0,2956,-1191,1309,917, $1084,0,-1559,-421,1212,84,778,360,887,0,678,218, $631,-109,-416,-173,178,-51,-211,0,327,-148,245, $58,-16,-34,-111,-126,-51,32,61,0,57,-2,-10,93, $-228,8,-51,-26,49,23,-98,-62,72,0,-54,-48,2, $-14,31,-12,-38,2,-18,28,19,-15,6,-22,11,0,8,8, $-4,-15,-9,6,2,-13,3,5,-8,16,6,-5,8,-18,8,0, $10,-20,1,14,-11,5,12,-3,1,-2,-2,8,2,10,0,-2, $-1,2,-3,0,-4,2,2,1,-5,2,-2,6,6,-4,4,0,0,-2, $1,4,2,0,0,-6/ data y1920/ $-31060,0,-2317,5845,-839,0,2959,-1259,1407,823, $1111,0,-1600,-445,1205,103,839,293,889,0,695,220, $616,-134,-424,-153,199,-57,-221,0,326,-122,236, $58,-23,-38,-119,-125,-62,43,61,0,55,0,-10,96, $-233,11,-46,-22,44,18,-101,-57,73,0,-54,-49,2, $-14,29,-13,-37,4,-16,28,19,-16,6,-22,11,0,7,8, $-3,-15,-9,6,2,-14,4,5,-7,17,6,-5,8,-19,8,0, $10,-20,1,14,-11,5,12,-3,1,-2,-2,9,2,10,0,-2, $-1,2,-3,0,-4,2,2,1,-5,2,-2,6,6,-4,4,0,0,-2, $1,4,3,0,0,-6/ data y1925/ $-30926,0,-2318,5817,-893,0,2969,-1334,1471,728, $1140,0,-1645,-462,1202,119,881,229,891,0,711,216, $601,-163,-426,-130,217,-70,-230,0,326,-96,226,58, $-28,-44,-125,-122,-69,51,61,0,54,3,-9,99,-238, $14,-40,-18,39,13,-103,-52,73,0,-54,-50,3,-14, $27,-14,-35,5,-14,29,19,-17,6,-21,11,0,7,8,-3, $-15,-9,6,2,-14,4,5,-7,17,7,-5,8,-19,8,0,10, $-20,1,14,-11,5,12,-3,1,-2,-2,9,2,10,0,-2,-1, $2,-3,0,-4,2,2,1,-5,2,-2,6,6,-4,4,0,0,-2,1, $4,3,0,0,-6/ data y1930/ $-30805,0,-2316,5808,-951,0,2980,-1424,1517,644, $1172,0,-1692,-480,1205,133,907,166,896,0,727,205, $584,-195,-422,-109,234,-90,-237,0,327,-72,218,60, $-32,-53,-131,-118,-74,58,60,0,53,4,-9,102,-242, $19,-32,-16,32,8,-104,-46,74,0,-54,-51,4,-15,25, $-14,-34,6,-12,29,18,-18,6,-20,11,0,7,8,-3,-15, $-9,5,2,-14,5,5,-6,18,8,-5,8,-19,8,0,10,-20, $1,14,-12,5,12,-3,1,-2,-2,9,3,10,0,-2,-2,2, $-3,0,-4,2,2,1,-5,2,-2,6,6,-4,4,0,0,-2,1,4, $3,0,0,-6/ data y1935/ $-30715,0,-2306,5812,-1018,0,2984,-1520,1550,586, $1206,0,-1740,-494,1215,146,918,101,903,0,744,188, $565,-226,-415,-90,249,-114,-241,0,329,-51,211,64, $-33,-64,-136,-115,-76,64,59,0,53,4,-8,104,-246, $25,-25,-15,25,4,-106,-40,74,0,-53,-52,4,-17,23, $-14,-33,7,-11,29,18,-19,6,-19,11,0,7,8,-3,-15, $-9,5,1,-15,6,5,-6,18,8,-5,7,-19,8,0,10,-20, $1,15,-12,5,11,-3,1,-3,-2,9,3,11,0,-2,-2,2, $-3,0,-4,2,2,1,-5,2,-2,6,6,-4,4,0,0,-1,2,4, $3,0,0,-6/ data y1940/ $-30654,0,-2292,5821,-1106,0,2981,-1614,1566,528, $1240,0,-1790,-499,1232,163,916,43,914,0,762,169, $550,-252,-405,-72,265,-141,-241,0,334,-33,208,71, $-33,-75,-141,-113,-76,69,57,0,54,4,-7,105,-249, $33,-18,-15,18,0,-107,-33,74,0,-53,-52,4,-18,20, $-14,-31,7,-9,29,17,-20,5,-19,11,0,7,8,-3,-14, $-10,5,1,-15,6,5,-5,19,9,-5,7,-19,8,0,10,-21, $1,15,-12,5,11,-3,1,-3,-2,9,3,11,1,-2,-2,2, $-3,0,-4,2,2,1,-5,2,-2,6,6,-4,4,0,0,-1,2,4, $3,0,0,-6/ data y1945/ $-30594,0,-2285,5810,-1244,0,2990,-1702,1578,477, $1282,0,-1834,-499,1255,186,913,-11,944,0,776,144, $544,-276,-421,-55,304,-178,-253,0,346,-12,194,95, $-20,-67,-142,-119,-82,82,59,0,57,6,6,100,-246, $16,-25,-9,21,-16,-104,-39,70,0,-40,-45,0,-18,0, $2,-29,6,-10,28,15,-17,29,-22,13,0,7,12,-8,-21, $-5,-12,9,-7,7,2,-10,18,7,3,2,-11,5,0,-21,-27, $1,17,-11,29,3,-9,16,4,-3,9,-4,6,-3,1,-4,8, $-3,0,11,5,1,1,2,-20,-5,-1,-1,-6,8,6,-1,-4, $-3,-2,5,0,-2,-2/ data y1950/ $-30554,0,-2250,5815,-1341,0,2998,-1810,1576,381, $1297,0,-1889,-476,1274,206,896,-46,954,0,792,136, $528,-278,-408,-37,303,-210,-240,0,349,3,211,103, $-20,-87,-147,-122,-76,80,54,0,57,-1,4,99,-247, $33,-16,-12,12,-12,-105,-30,65,0,-55,-35,2,-17, $1,0,-40,10,-7,36,5,-18,19,-16,22,0,15,5,-4, $-22,-1,0,11,-21,15,-8,-13,17,5,-4,-1,-17,3,0, $-7,-24,-1,19,-25,12,10,2,5,2,-5,8,-2,8,3,-11, $8,-7,-8,0,4,13,-1,-2,13,-10,-4,2,4,-3,12,6, $3,-3,2,6,10,11,3,8/ data y1955/ $-30500,0,-2215,5820,-1440,0,3003,-1898,1581,291, $1302,0,-1944,-462,1288,216,882,-83,958,0,796,133, $510,-274,-397,-23,290,-230,-229,0,360,15,230,110, $-23,-98,-152,-121,-69,78,47,0,57,-9,3,96,-247, $48,-8,-16,7,-12,-107,-24,65,0,-56,-50,2,-24,10, $-4,-32,8,-11,28,9,-20,18,-18,11,0,9,10,-6,-15, $-14,5,6,-23,10,3,-7,23,6,-4,9,-13,4,0,9,-11, $-4,12,-5,7,2,6,4,-2,1,10,2,7,2,-6,5,5,-3,0, $-5,-4,-1,0,2,-8,-3,-2,7,-4,4,1,-2,-3,6,7,-2, $-1,0,-3/ data y1960/ $-30421,0,-2169,5791,-1555,0,3002,-1967,1590,206, $1302,0,-1992,-414,1289,224,878,-130,957,0,800, $135,504,-278,-394,3,269,-255,-222,0,362,16,242, $125,-26,-117,-156,-114,-63,81,46,0,58,-10,1,99, $-237,60,-1,-20,-2,-11,-113,-17,67,0,-56,-55,5, $-28,15,-6,-32,7,-7,23,17,-18,8,-17,15,0,6,11, $-4,-14,-11,7,2,-18,10,4,-5,23,10,1,8,-20,4,0, $6,-18,0,12,-9,2,1,0,4,-3,-1,9,-2,8,3,0,-1, $5,1,0,-3,4,4,1,0,0,-1,2,4,-5,6,1,1,-1,-1, $6,2,0,0,-7/ data y1965/ $-30334,0,-2119,5776,-1662,0,2997,-2016,1594,114, $1297,0,-2038,-404,1292,240,856,-165,957,0,804, $148,479,-269,-390,13,252,-269,-219,0,358,19,254, $128,-31,-126,-157,-97,-62,81,45,0,61,-11,8,100, $-228,68,4,-32,1,-8,-111,-7,75,0,-57,-61,4,-27, $13,-2,-26,6,-6,26,13,-23,1,-12,13,0,5,7,-4, $-12,-14,9,0,-16,8,4,-1,24,11,-3,4,-17,8,0,10, $-22,2,15,-13,7,10,-4,-1,-5,-1,10,5,10,1,-4, $-2,1,-2,0,-3,2,2,1,-5,2,-2,6,4,-4,4,0,0,-2, $2,3,2,0,0,-6/ data y1970/ $-30220,0,-2068,5737,-1781,0,3000,-2047,1611,25, $1287,0,-2091,-366,1278,251,838,-196,952,0,800, $167,461,-266,-395,26,234,-279,-216,0,359,26,262, $139,-42,-139,-160,-91,-56,83,43,0,64,-12,15,100, $-212,72,2,-37,3,-6,-112,1,72,0,-57,-70,1,-27, $14,-4,-22,8,-2,23,13,-23,-2,-11,14,0,6,7,-2, $-15,-13,6,-3,-17,5,6,0,21,11,-6,3,-16,8,0,10, $-21,2,16,-12,6,10,-4,-1,-5,0,10,3,11,1,-2,-1, $1,-3,0,-3,1,2,1,-5,3,-1,4,6,-4,4,0,1,-1,0, $3,3,1,-1,-4/ data y1975/ $-30100,0,-2013,5675,-1902,0,3010,-2067,1632,-68, $1276,0,-2144,-333,1260,262,830,-223,946,0,791, $191,438,-265,-405,39,216,-288,-218,0,356,31,264, $148,-59,-152,-159,-83,-49,88,45,0,66,-13,28,99, $-198,75,1,-41,6,-4,-111,11,71,0,-56,-77,1,-26, $16,-5,-14,10,0,22,12,-23,-5,-12,14,0,6,6,-1, $-16,-12,4,-8,-19,4,6,0,18,10,-10,1,-17,7,0, $10,-21,2,16,-12,7,10,-4,-1,-5,-1,10,4,11,1, $-3,-2,1,-3,0,-3,1,2,1,-5,3,-2,4,5,-4,4,-1, $1,-1,0,3,3,1,-1,-5/ data y1980/ $-29992,0,-1956,5604,-1997,0,3027,-2129,1663,-200, $1281,0,-2180,-336,1251,271,833,-252,938,0,782, $212,398,-257,-419,53,199,-297,-218,0,357,46,261, $150,-74,-151,-162,-78,-48,92,48,0,66,-15,42,93, $-192,71,4,-43,14,-2,-108,17,72,0,-59,-82,2,-27, $21,-5,-12,16,1,18,11,-23,-2,-10,18,0,6,7,0, $-18,-11,4,-7,-22,4,9,3,16,6,-13,-1,-15,5,0, $10,-21,1,16,-12,9,9,-5,-3,-6,-1,9,7,10,2,-6, $-5,2,-4,0,-4,1,2,0,-5,3,-2,6,5,-4,3,0,1,-1, $2,4,3,0,0,-6/ data y1985/ $-29873,0,-1905,5500,-2072,0,3044,-2197,1687,-306, $1296,0,-2208,-310,1247,284,829,-297,936,0,780, $232,361,-249,-424,69,170,-297,-214,0,355,47,253, $150,-93,-154,-164,-75,-46,95,53,0,65,-16,51,88, $-185,69,4,-48,16,-1,-102,21,74,0,-62,-83,3,-27, $24,-2,-6,20,4,17,10,-23,0,-7,21,0,6,8,0,-19, $-11,5,-9,-23,4,11,4,14,4,-15,-4,-11,5,0,10, $-21,1,15,-12,9,9,-6,-3,-6,-1,9,7,9,1,-7,-5, $2,-4,0,-4,1,3,0,-5,3,-2,6,5,-4,3,0,1,-1,2, $4,3,0,0,-6/ data y1990/ $-29775,0,-1848,5406,-2131,0,3059,-2279,1686,-373, $1314,0,-2239,-284,1248,293,802,-352,939,0,780, 247,325,-240,-423,84,141,-299,-214,0,353,46,245, $154,-109,-153,-165,-69,-36,97,61,0,65,-16,59,82, $-178,69,3,-52,18,1,-96,24,77,0,-64,-80,2,-26, $26,0,-1,21,5,17,9,-23,0,-4,23,0,5,10,-1,-19, $-10,6,-12,-22,3,12,4,12,2,-16,-6,-10,4,0,9, $-20,1,15,-12,11,9,-7,-4,-7,-2,9,7,8,1,-7,-6, $2,-3,0,-4,2,2,1,-5,3,-2,6,4,-4,3,0,1,-2,3, $3,3,-1,0,-6/ data y1995/ $-29692,0,-1784,5306,-2200,0,3070,-2366,1681,-413, $1335,0,-2267,-262,1249,302,759,-427,940,0,780, $262,290,-236,-418,97,122,-306,-214,0,352,46,235, $165,-118,-143,-166,-55,-17,107,68,0,67,-17,68, $72,-170,67,-1,-58,19,1,-93,36,77,0,-72,-69,1, $-25,28,4,5,24,4,17,8,-24,-2,-6,25,0,6,11,-6, $-21,-9,8,-14,-23,9,15,6,11,-5,-16,-7,-4,4,0, $9,-20,3,15,-10,12,8,-6,-8,-8,-1,8,10,5,-2,-8, $-8,3,-3,0,-6,1,2,0,-4,4,-1,5,4,-5,2,-1,2, $-2,5,1,1,-2,0,-7/ data y2000/ $-29619.4,0,-1728.2,5186.1,-2267.7,0,3068.4,-2481.6, $1670.9,-458.0,1339.6,0,-2288.0,-227.6,1252.1,293.4, $714.5,-491.1,932.3,0,786.8,272.6,250.0,-231.9, $-403.0,119.8,111.3,-303.8,-218.8,0,351.4,43.8,222.3, $171.9,-130.4,-133.1,-168.6,-39.3,-12.9,106.3,72.3,0, $68.2,-17.4,74.2,63.7,-160.9,65.1,-5.9,-61.2,16.9, $0.7,-90.4,43.8,79.0,0,-74.0,-64.6,0.0,-24.2,33.3, $6.2,9.1,24.0,6.9,14.8,7.3,-25.4,-1.2,-5.8,24.4,0, $6.6,11.9,-9.2,-21.5,-7.9,8.5,-16.6,-21.5,9.1,15.5, $7.0,8.9,-7.9,-14.9,-7.0,-2.1,5.0,0,9.4,-19.7,3.0, $13.4,-8.4,12.5,6.3,-6.2,-8.9,-8.4,-1.5,8.4,9.3, $3.8,-4.3,-8.2,-8.2,4.8,-2.6,0,-6.0,1.7,1.7,0.0, $-3.1,4.0,-0.5,4.9,3.7,-5.9,1.0,-1.2,2.0,-2.9,4.2, $0.2,0.3,-2.2,-1.1,-7.4/ data y2005/ $-29556.8,0,-1671.8,5080.0,-2340.5,0,3047.0,-2594.9, $1656.9,-516.7,1335.7,0,-2305.3,-200.4,1246.8,269.3, $674.4,-524.5,919.8,0,798.2,281.4,211.5,-225.8, $-379.5,145.7,100.2,-304.7,-227.6,0,354.4,42.7,208.8, $179.8,-136.6,-123.0,-168.3,-19.5,-14.1,103.6,72.9,0, $69.6,-20.2,76.6,54.7,-151.1,63.7,-15.0,-63.4,14.7, $0.0,-86.4,50.3,79.8,0,-74.4,-61.4,-1.4,-22.5,38.6, $6.9,12.3,25.4,9.4,10.9,5.5,-26.4,2.0,-4.8,24.8,0, $7.7,11.2,-11.4,-21.0,-6.8,9.7,-18.0,-19.8,10.0, $16.1,9.4,7.7,-11.4,-12.8,-5.0,-0.1,5.6,0,9.8, $-20.1,3.6,12.9,-7.0,12.7,5.0,-6.7,-10.8,-8.1,-1.3, $8.1,8.7,2.9,-6.7,-7.9,-9.2,5.9,-2.2,0,-6.3,2.4, $1.6,0.2,-2.5,4.4,-0.1,4.7,3.0,-6.5,0.3,-1.0,2.1, $-3.4,3.9,-0.9,-0.1,-2.3,-2.2,-8.0, $2.7,0,-1.7,0.1,-1.9,1.3,1.5,-0.9,-0.1,-2.6,0.1, $0.9,-0.7,-0.7,0.7,-2.8,1.7,-0.9,0.1,-1.2,1.2,-1.9, $4.0,-0.9,-2.2,0,-0.3,-0.4,0.2,0.3,0.9,2.5,-0.2, $-2.6,0.9,0.7,-0.5,0.3,0.3,0.0,-0.3,0.0,-0.4,0.3, $-0.1,-0.9,-0.2,-0.4,-0.4,0.8,-0.2,0,-0.9,-0.9,0.3, $0.2,0.1,1.8,-0.4,-0.4,1.3,-1.0,-0.4,-0.1,0.7,0.7, $-0.4,0.3,0.3,0.6,-0.1,0.3,0.4,-0.2,0.0,-0.5,0.1,-0.9/ end c______________________________________________________________ c======================================================================= c UNIT K: DECODING ROUTINES FOR COMMAND KIT c======================================================================= subroutine getarr(code, values, nwant, nfound) c$$$$ calls nothing c c Extracts an array of numbers from the input store. That store is c a large array in common /store/ which has been filled earlier c c code A 4-bye identifying code. Only lines in the input store c beginning with this code are scanned for numbers. c values the real output array of values found. c nwant the maximum number of numbers expected . c nfound the number of numbers actually found in the input. c If the line contains fewer than nwant values, this is the c value returned in nfound. If an error is discovered c nfound=-n, where n is the number of numbers properly c decoded. If there are no numbers after the codeword c nfound=0. Finally, if the code is absent from the store c nfound=-99 and the array is left undisturbed. c parameter (inmx=200) character *80 input(inmx),line common /store/ input common /ndict/ iecho,nin,memory,kread(inmx) dimension values(*) character local*40, code*4, char*80 c c Read the store in normal order ce do 1010 lin=nin, 1, -1 line=input(lin) c Check code and if line is fresh or old if (code .eq. line(1:4)) then if (iecho .ge. 1) write(*,'(2a)')'==> ',line c Increment read count; copy tail into char, discarding stuff after % kread(lin)=kread(lin) + 1 n1=index(line, ' ')+1 n2=index(line, '%')-1 if (n2 .le. 0) n2=80 char=line(n1:n2) kn=n2 - n1 + 1 goto 1020 endif 1010 continue c Code word not found nfound=-99 return c 1020 continue k1=1 c Up to nwant numbers are sought do 1800 nval=1, nwant do 1100 k=k1, kn if (char(k:k) .ne. ' ') goto 1200 1100 continue nfound=nval-1 return 1200 do 1300 l=k, kn if (char(l:l).eq. ',' .or. char(l:l).eq.' ') goto 1500 1300 continue 1500 if (index(' Ee+-.0123456789', char(k:k)) .eq. 0) goto 2000 if (k-l+1 .gt. 40) goto 2010 k2=l+1 local=char(k:l-1) read (local, '(f40.0)', err=2000) values(nval) 1800 k1=k2 nval=1 - nwant 1900 nfound=1 - nval return c 2000 print '(a)',' ', $'>>> Unreadable numbers in this input line:',line goto 1900 2010 print '(a)',' ', $'>>> Program can not read a number of more than 40 characters:', $char(k:l-1) goto 1900 end c______________________________________________________________ subroutine getone(code, value, nfound) c$$$$ calls getarr c c Extracts a single number from the input store. That store is c a large array in common /store/ which has been filled earlier c c code A 4-bye identifying code. Only lines in the input store c beginning with this code are scanned for numbers. c value the real output variable containing the desired number. c nfound is 1 if a number is successfully read in; it is zero c the number is absent or unreadable. nfound = -99 if the c code is absent from the input store. c character*4 code dimension v(1) ce call getarr(code, v, 1, nfound) if (nfound .eq. 1) value=v(1) return end c______________________________________________________________ subroutine getchr(code, char, nfound) c$$$$ calls nothing c Extracts a single character variable from the input store. That c store is a large array in common /store/ which has been filled c earlier c c code A 4-bye identifying code. Only lines in the input store c beginning with this code are scanned. c char the character output variable containing the desired string. c nfound is 1 if a string is successfully read in; it is zero c the line is blank after the code. nfound = -99 if the c code is absent from the input store. c parameter (inmx=200) character *80 input(inmx),line, char*(*), code*4 common /store/ input common /ndict/ iecho,nin,memory,kread(inmx) c Inspect the store in normal order ce do 1010 lin=nin, 1, -1 line=input(lin) c Check code and if line is fresh or old if (code .eq. line(1:4)) then if (iecho .ge. 1) write(*,'(2a)')'==> ',line c Increment read count kread(lin)=kread(lin) + 1 goto 1020 endif 1010 continue c Code word not found nfound=-99 return c 1020 continue c n1=index(line, ' ')+1 n2=index(line, '%')-1 if (n2 .le. 0) n2=80 c Save in char everything after 1st blank to just before % do 1200 k=n1, n2 if (line(k:k) .ne. ' ') then char=line(k:n2) nfound=1 return endif 1200 continue c c Blank field after code word nfound=0 return end c______________________________________________________________ subroutine clear(code) c$$$$ calls nothing c c Removes all command entries identified by code. c code A 4-bye identifying code. Only lines in the input store c beginning with this code are scanned. c if byte 5 is '1' clear only one command not all c parameter (inmx=200) character *80 input(inmx),line, code*5 common /store/ input common /ndict/ iecho,nin,memory,kread(inmx) c c Read the store in normal order ce do 1010 lin= nin, 1, -1 line=input(lin) c Check code and clear the line if present through stack if (code(1:4) .eq. line(1:4)) then input(lin)=' cleared' if (code(5:5).eq.'1') return endif 1010 continue return end c______________________________________________________________ c======================================================================= c UNIT Y: GENERATE AND SUM Ylm c======================================================================= subroutine shsum(ra, lmx,blm, br,bth,bla,V) c$$$$ calls allm c Evaluates a vector-valued function -grad V, and V, where V is a c harmonic function with internal (lmx>0) or external (lmx<0) c sources and has complex spherical harmonic coefficients blm, c fully normalized. Thus c when lmx>0: V = sum{l,m} [b_lm*(1/r)**(l+1) * Ylm] c when lmx<0: V = sum{l,m} [b_lm*r**l * Ylm] c The coefficients are stored serially c c ra 3-vector defining obs point c lmx spherical harmonic degree of series. c if lmx < 0 fields are of external origin, and l(max)=-lmx! c blm complex coefficients in a fully normalized spherical harmonic c series for V. Stored serially. c values are given only for 0 .le. m .le. lmax. c br returned radial component of -grad V. c bth returned value of theta (colatitude) component of -grad V c bla returned value of lambda (longitude) component of -grad V c V returned value of scalar potential V c parameter (mxdeg=96) c dimension r(3),ra(3) complex dx,slm,ex(mxdeg), blm(*) double precision plm,dlm,cosco dimension plm(mxdeg+1),dlm(mxdeg+1) c data near/0/ c ce lmax=abs(lmx) if (lmax .ge. mxdeg) then print '(a,i5/a)','>>> lmax too large for shsum:', lmax, $ '>>> Reduce lmax or recompile shsum with bigger mxdeg' stop endif c Form unit vector rr=sqrt(ra(1)**2 + ra(2)**2 + ra(3)**2) r(1)=ra(1)/rr r(2)=ra(2)/rr r(3)=ra(3)/rr c sinth=sqrt(r(1)**2 + r(2)**2) if (sinth .eq. 0.0) then r(1)=2.0e-6 sinth=r(1) if (near .lt. 1) print '(a)', $ '>>> One or more observer sites exactly at pole:', $ '>>> Magmap moves such sites 7 meters' near=near + 1 elseif (sinth .lt. 2.0e-6) then r(1)=r(1)*2.0e-6/sinth r(2)=r(2)*2.0e-6/sinth sinth=2.0e-6 endif c cosco=r(3) if (abs(r(3)).gt. 0.99) $ cosco=sqrt(1.0d00 - sinth**2)*sign(1.0, r(3)) c c Calculate and save exp(i*m*phi) ex(1)=1.0 dx=cmplx(r(1), r(2))/sinth do 1300 m=1, lmax ex(m+1)=dx*ex(m) 1300 continue c c Sum the series for -grad V and V br =0.0 bth=0.0 bla=0.0 V=0.0 j=0 do 1700 l=0, lmax call allm(l, cosco, plm,dlm) c Chooses external/internal part by sign of lmx n=l + 2 if (lmx.lt. 0) n=1 - l rpow= rr**n do 1600 m=0, l j=j + 1 slm=ex(m+1)*blm(j) / rpow em=m V =V + rr* plm(m+1)*real(slm) br =br + (n - 1.0) * plm(m+1)*real(slm) bth=bth - dlm(m+1)*real(slm) bla=bla - plm(m+1)*real(cmplx(0.0, em)*slm)/sinth 1600 continue 1700 continue c return end c_______________________________________________________________________ subroutine allm(l, amu, plm,dlm) c$$$$ calls nothing c For degree l finds fully-normalized associated Legendre functions c of order m= 0, 1, 2, ... l for the argument amu=cos(theta) and c corresponding theta derivatives (theta=colatitude). c Results of order m held in plm(m+1), dlm(m+1) c c plm(m+1)=P{l,m}(x)=t(m)*(-1)**m*sqrt((2*l+1)*(l-m)!/(l+m)!)*Plm(x) c where c Plm(x)=(1/2**l*l!)*(1-x**2)**(m/2)*(d/dx)**(l+m)(x**2-1)**l c and t(0)=1/sqrt(4*pi), t(m)=1/sqrt(2*pi), m>0. c c Note that 0 =< l =< 100, unless parameter lmx is changed c implicit double precision (a-h, o-z) parameter (lmx=100, nmx=2*lmx, twopi=2*3.14159265358979324d0) dimension plm(l+1),dlm(l+1),rootn(0:nmx),tlm(lmx+1) save init,rootn data init/0/ c c Fill array with square-roots of the integers ce if (init .eq. 0) then do 1000 n=0, nmx en=n rootn(n)=sqrt(en) 1000 continue init=1 endif c c Check degree constraints if (l .gt. lmx) then print '(a,i5)','>>> Degree too great for ALLM: l >',lmx stop elseif (l .le. 0) then plm(1)=1.0/sqrt(2.0*twopi) return endif c c Produce normalization constant c=(l + 0.5d0)/twopi twol=2.0*l do 1100 x=2.0, twol, 2.0 c=c - c/x 1100 continue s=sqrt(1.0 - amu**2) twocot =-2.0*amu/s tlm(l+1)= 2.0*sqrt(c)*(-s)**l tlm(l) = tlm(l+1)*twocot*l/sqrt(twol) c Recur downwards to m=0 do 1500 m=l-2, 0, -1 tlm(m+1)=((m+1)*twocot*tlm(m+2) - $ rootn(l+m+2)*rootn(l-m-1)*tlm(m+3))/(rootn(l+m+1)*rootn(l-m)) 1500 continue c c Calculate the derivative wrt theta cot=amu/s do 1600 m=0, l-1 plm(m+1)=tlm(m+1) dlm(m+1)=rootn(l+m+1)*rootn(l-m)*tlm(m+2)+m*cot*tlm(m+1) 1600 continue dlm(l+1)=l*cot*tlm(l+1) plm(l+1)=tlm(l+1) plm(1)=tlm(1)/2.0 dlm(1)=dlm(1)/2.0 return end c_______________________________________________________________________ c======================================================================= c UNIT P: SINGLE POINT EVALUATION c======================================================================= subroutine point c$$$$ calls getarr shsum clear c Calculates a field model at a specified geographical point parameter (mxdeg=96, mxb=(mxdeg+1)*(mxdeg+2)/2, mxa=mxb/10) complex blm,alm common /coeff/ lmax,lext,blm(mxb),alm(mxa) common /io/ inp,iout,iprint common /earth/ rad(4) dimension site(3),r(3) character*64 loci c data radian/0.01745329/ c ce a=1000.0*rad(1) b=1000.0*rad(2) if (rad(4) .le. 0.7) then a=1000.0*rad(3) b=a endif if (a .eq. b) then if (iprint.ge.0) print '(a)',' ',' NOTE:', $ ' Ellipticity is always omitted for core evaluations', $ ' Hence geocentric coordinate system is in use',' ' if (iprint.ge.0) $ print '(a,f9.1)',' Standard evaluation radius (km) ',0.001*a endif c call getchr('site', loci, nsite) if (nsite .le. 0) return open(unit=12, file=loci, status='OLD',err=1100) if (iprint.ge.0) $print '(a)',' ',' Evaluation sites read from:',loci nsite=99 goto 1200 1100 continue 1200 do 1800 k=1, 5500 if (nsite .ne. 99) then call getarr('site', site, 3, nsite) c Clear just one site command call clear('site1') height=0 if (nsite .eq. 3) height=site(3) else read (12,*, end=1810, err=2000) site(1),site(2),site(3) height=site(3) endif if (nsite .le. 0) return if (k .eq. 1 .and.iprint.ge.0) print '(a)', ' ', $ ' Site values automatically sent to file fort.9. ', $ ' Each line in the file consists of:', $ ' lat, long, height, X, Y, Z, F, D, I',' ' alat=site(1) along=site(2) c Coordinate transformations from Langel, pp264-267 in Jacobs c "GEOMAGNETISM", 1987 c c Convert geodetic coordinates to geocentric ones sinla=sin(radian*alat) cosla=cos(radian*alat) r0=sqrt((a*cosla)**2 + (b*sinla)**2) x1=(a**2/r0 + height)*cosla y1=(b**2/r0 + height)*sinla r1=sqrt(x1**2 + y1**2) rabs=r1/a cost=y1/r1 sint=x1/r1 c c Find field elements in geocentric frame c Recall rad(3) is the reference radius scale=rad(1)*rabs / rad(3) r(3)=scale*cost r(1)=scale*sint*cos(radian*along) r(2)=scale*sint*sin(radian*along) call shsum(r, lmax, blm, Br,Bt,Bl,V) if (lext.gt. 0) then call shsum(r, -lext, alm, exBr,exBt,exBl,exV) Br=Br + exBr Bt=Bt + exBt Bl=Bl + exBl V =V + exV endif c c Convert field elements to geodetic frame sins=sinla*sint - cos(radian*alat)*cost coss=sqrt(1.0 - sins**2) F=sqrt(Br**2 + Bt**2 + Bl**2) X=-Bt*coss - Br*sins Y=Bl Z=Bt*sins - Br*coss decl=atan2(Y, X+1e-18)/radian ainc=atan2(Z, sqrt(X**2 + Y**2))/radian c c Print it out if (F .lt. 100000.0) then if (iprint.ge.0) $ print '(/(a18,f10.2,a6,f10.2,a6,f10.2,a))', $' Site: lat ',alat,' long',along,' alt',height,' m', $' Elements (nT): X', X,' Y', Y,' Z',Z,' ', $' F', F,' Decl', decl,' Inc',ainc,' ' else if (iprint.ge.0) $ print '(/ a18,f12.2,a6,f12.2,a6,f12.2,a/ $ (+1p,a18,g12.4,a6,g12.4,a6,g12.4,a))', $' Site: lat ',alat,' long',along,' alt',height,' m', $' Elements (nT): X', X,' Y', Y,' Z',Z,' ', $' F', F,' Decl', decl,' Inc',ainc,' ' endif write (8,'(f6.2,f8.2,1p,e11.4,0p,4f9.1,2f8.2)') $ alat,along,height,X,Y,Z,F,decl,ainc 1800 continue print '(a)','>>> Some sites omitted from evaluation' 1810 close(unit=12) return c c Error out 2000 print '(a,i7)','>>> Error reading site file at line:',k close (unit=12) return end c_______________________________________________________________________ subroutine plot(m1,n1, iproj, name) c$$$$ calls launch getchr projg c Writes a contour or color command file and executes it c Also writes the coast outline file in mapped coordinates c m1,n1 array size; iproj projection type; name filename of array save kgrph character*(*) name character*64 shore, oldsho, script, tscrip, costid*8 character*64 syline, syfile*10 dimension geo(2), u(2) common /io/ inp,iout,iprint common /scan/ m,n,umin,vmin,du,factor,hash common /world/ c(3217),nworld parameter (rd=.0174533, deg=rd) data bdum/1.0e8/, u1,u2/2*1e9/ data script/' '/,nshor/0/, oldhsh/0/,idpts/0/ data oldsho/' '/ c c Run plotfile in script - final act if (m1 .eq. 0) then call launch(kgrph, script) return endif c c Decide if a plot is called for call getchr('color', tscrip, nscr1) call getchr('plot', tscrip, nscr2) c Color map requested (kgrph=1) nscr=max(nscr1, min(0,nscr2)) lk=kgrph if (nscr .ge. 0) kgrph =1 if (nscr .lt. 0) then call getchr('contour', tscrip, nscr) c Contour map requested (kgrph=0) if (nscr .ge. 0) kgrph=0 endif c No plot required if (nscr .lt. 0) return c if (nscr .eq. 0) tscrip='tmp.plt' c New plotfile requested if (tscrip .ne. script) then c Run preceeding plotfile then open a new plotfile call launch(lk, script) script=tscrip open(unit=10, name=script) endif c c Create a file of mapped coastline to be superimposed on contor map call getchr('coast', shore, nsh) if (nsh .ge. 0) then c Don't create a new coast file if projection is unchanged if (hash.eq.oldhsh .and. shore.eq.oldsho) goto 1200 c if (nsh .gt. 0) open(unit=8, file=shore, status='old',err=2200) nshor=nshor+1 oldsho=shore write(costid,'(a,i2)') 'coast',nshor+9 open (unit=11, file=costid) c crit=0.2 if (iproj .eq. 1) crit=0.2/deg ncoast=0 do 1100 j=1, nworld c Read external file or unpack internal coastline if (nsh .eq. 1) then read (8,*, err=2000, end=1110) geo else geo(1)=0.1*c(j) geo(2)=1000.0*mod(abs(c(j)), 1.0) endif c Transform lat-long data to x-y plane call projg(iproj, geo, u) dist=abs(u1-u(1)) + abs(u2-u(2)) c If jump is too large insert a pen-up command if (dist .gt. crit .and. $ geo(1).lt.90.0 .and. u1.lt. 90.0) then ncoast=ncoast + 1 write (11,'(a)') '99 370' endif ncoast=ncoast + 1 write (11,'(2f9.4)') u u1=u(1) u2=u(2) 1100 continue c 1110 if (nsh .eq. 1) close (unit=8) close (unit=11) if (iprint.ge.0) print '(a,i5,a)', $ ' A coastline has been added to the map of',ncoast,' points' endif 1200 continue c c Create files of mapped sample points to be mapped and superimposed do 1400 nfile=1, 100 call getchr('symbols', syline, mysy) call clear( 'symbols') c Delete symbol file if blank field or when projection changes if (mysy.eq.0 .or. (oldhsh*mysy.lt.0 .and. hash.ne.oldhsh)) then print '(/a)',' Symbol data removed from this plot' print '(a)',' To retain symbol data in a new projection, '// $ 'repeat symbol command(s)' write(10,'(a)')'symbols' goto 1410 elseif (mysy .lt. 0) then goto 1410 elseif (mysy .gt. 0) then c Create new symbol file mapped to new coordinates isy=index(syline, ' ') close (unit=8) open (unit=8, file=syline(1:isy), status='OLD',err=1390) idpts=idpts+1 write(syfile,'(a,i2)') 'sym',idpts+9 close (unit=11) open (unit=11, file=syfile) print '(2a)',' Transforming coords in file: ',syline(1:isy) c c Read from symbol file, map point and write it out to the new file do 1300 j=1, nworld read (8,*, err=2000, end=1310) geo call projg(iproj, geo, u) write (11,'(2f9.4)') u 1300 continue c 1310 close (unit=8) close (unit=11) write(10,'(3a)') 'symbols ',syfile,syline(isy : 64) endif goto 1400 c 1390 print '(a,a)','>>> Unable to open file: ',syline(1:isy) 1400 continue c c Save id of current map projection 1410 oldhsh=hash c c Create a file to draw the map and execute it: kgrph=0 means c a contour map, =1 means color goto (1500, 1600), kgrph+1 c c Write command file for contour c Draw a bounding circle or ellipse on the map 1500 if (iproj .gt. 1 .and. iproj.ne.5) then write (10,'(a)') 'axes -1 1 -1 1','file *','read 181' write (10,'(10f8.4)') $ (cos(rd*j)/factor, sin(rd*j)/factor,j=0,360,2) endif write (10,'(a,a/a,3i5)') 'file ',name,'read ',m,n,1 write (10,'(a,2f9.1)') 'null ',1-bdum/1000, bdum/1000-1 write(10,'(a,4f9.4,a)')'axes ',umin,-umin,vmin,-vmin if (nsh .ge. 0) write (10,'(a/a/a,i6)') $'Color 9','file '//costid, 'read ',ncoast write (10,'(a)') 'axes 0 0 0 0', 'Color 0','letter 0.06' write (10,'(a)') 'width 5','fancy 0','plot ' return c c This is the color map code c Draw a bounding circle or ellipse on the map 1600 if (iproj .gt. 1 .and. iproj.ne.5) $ write (10,'(a)') 'border mask' write (10,'(a,a/a,3i5)') 'file ',name,'read ',m,n,1 write(10,'(a,4f9.4,a)')'axes ',umin,-umin,vmin,-vmin,' hidden' if (nsh .ge. 0) write (10,'(a,a)') 'lines ',costid write (10,'(a)') 'width 5','outline 0','palette 2', $'plot ' return c c Error return 2000 print '(2a/a,i6)','>>> Unreadable data in file',shore, $ ' on line ',j return 2200 print '(a,a)','>>> Unable to open file:',shore goto 1200 c end c_______________________________________________________________________