cd /home/SHARE manxtE.txt prog/prtps -p lib/man/manxyE.txt | lpr ~~~~~~~~~~ Conversion for Map Projection ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ NAME xyconv - Lat.Lon. to Japanese Transverse Mercator (Old Tokyo Datum) X,Y nxycnv - Lat.Lon. to Japanese Transverse Mercator (New JGD2000 Datum) X,Y ikconv - Japanese Transverse Mercator (Old Tokyo Datum) X,Y to Lat.Lon. nikcnv - Japanese Transverse Mercator (New JGD2000 Datum) X,Y to Lat.Lon. utm - Lat.Lon. to Universal Transverse Mercator X,Y utmik - Universal Transverse Mercator X,Y to Lat.Lon. [i: Ido = Latitude ; k: Keido = Longitude] (Note) For UTM, X¡á0. on the equator, Y¡á500.(km) on the central meridian. SYNOPSIS call xyconv(gi, gk, fi, fk, x, y) call nxycnv(gi, gk, fi, fk, x, y) call utm(nc, fi, fk, x, y) gi, gk Lat. and Lon. (in minutes) of origin nc Zone number of UTM coordinaets (1-60) for Old Tokyo Datum, or plus 200 (201-260) for New JGD2000 Datum fi, fk Lat. and Lon. (in minutes) of the point to be calculated x, y Calculated results. Northword(X) and Eastword(Y) coordinates (in km) values are given. call ikconv(gi, gk, x, y, fi, fk) call nikcnv(gi, gk, x, y, fi, fk) call utmik(nc, x, y, fi, fk) gi, gk Lat. and Lon. (in minutes) of origin nc Zone number of UTM coordinaets (1-60) for Old Tokyo Datum, or plus 200 (201-260) for New JGD2000 Datum x, y Northword(X) and Eastword(Y) coordinates (in km) of the point to be calculated. fi, fk Calculated results. Lat. and Lon. (in minutes) are given. NAME cvinit / cvdinit - Specify projection method cviken / cvdiken - Convert Lat.Lon. into plane rectangular coordinates cvenik / cvdenik - Convert plane rectangular coordinates into Lat.Lon. Two models of Earth's ellipsoid are implemented: 1. Bessel's ellipsoid : a = 6377.397155 km , flattening = 1/299.152813 2. GRS80 ellipsoid : a = 6378.137 km , flattening = 1/298.257222 Following projection numbers are available for Bessel's ellipsoid, while the numbers for GRS80 ellipsoid are plus 200 for each projection. 0 Japanese Transverse Mercator [JTM] 1-60 Universal Transverse Mercator [UTM] (with zone number) 61 Universal Polar Stereographic Projection [UPS] North Pole 62 Universal Polar Stereographic Projection [UPS] South Pole 65 Same as UTM but non-standard Central Meridian 70 Mercator projection 71 Lambert Conformal Conic projection with 1 Standard Pararrel 72 Lambert Conformal Conic projection with 2 Standard Pararrels 100 Lambert Azimuthal Equal Area projection from a sphere with the same surface area as the ellipsoid 109 Lambert Azimuthal Equal Area projection from a sphere with the same eqatorial radius as the ellipsoid (Note) Although X¡áY¡á0. at the origin in general, for UTM, X=0. on the equator, and Y=500.(km) on the central meridian, for UPS, X¡áY¡á2000.(km) at the pole. SYNOPSIS call cvinit(nc, gi, gk, sp1, sp2) call cvdinit(nc, dgi, dgk, dp1, dp2) nc Projection number gi, gk Lat. and Lon. (in minutes) of origin (X¡áY¡á0.) (not refered in cases of UTM/UPS) dgi, dgk [ditto] sp1, sp2 Lat. (in minutes) of two standard pararells (not refered in cases other than projection 72,272) dp1, dp2 [ditto] call cviken(fi, fk, ye, xn) call cvdiken(dfi, dfk, dye, dxn) fi, fk Lat. and Lon. (in minutes) of the point to be calculated dfi, dfk [ditto] ye, xn Calculated results. Eastword(Y) and Northword(X) coordinates (in km) values are given. dye, dxn [ditto] call cvenik(ye, xn, fi, fk) call cvdenik(dye, dxn, dfi, dfk) ye, xn Eastword(Y) and Northword(X) coordinates (in km) of the point to be calculated. dye, dxn [ditto] fi, fk Calculated results. Lat. and Lon. (in minutes) are given. dfi, dfk [ditto] (Note) Be aware that the sequence of coordinates values (X, Y) arguments in cviken / cvenik / cvdiken / cvdenik is different from those in utm / utmik / xyconv / ikconv / nxycnv / nikcnv . -------------------------------- Geodetic Coordinates Conversion ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ NAME xw84t / xw84td - From World Geodetic System to Tokyo Datum xtw84 / xtw84d - From Tokyo Datum to World Geodetic System after Kanazawa(1988): "Calculation formula for the relationship between WGS-84 and Tokyo Datum indicated on Japanese charts. Data Rept. Hydrogr. Obs., Ser. Satellite Geodesy, no.1, 76-81". SYNOPSIS call xw84t (wlat, wlon, walt, tlat, tlon, talt) call xw84td(dwlat, dwlon, dwalt, dtlat, dtlon, dtalt) wlat, wlon WGS Lat. and Lon. (in minutes) to be converted dwlat,dwlon [ditto] walt WGS altitude (in m) to be converted dwalt [ditto] tlat, tlon Calculated results. Tokyo Datum Lat. and Lon. (in minutes) are given. dtlat,dtlon [ditto] talt Calculated result. Tokyo Datum altitude (in m) is given. dtalt [ditto] call xtw84 (tlat, tlon, talt, wlat, wlon, walt) call xtw84d(dtlat, dtlon, dtalt, dwlat, dwlon, dwalt) tlat, tlon Tokyo Datum Lat. and Lon. (in minutes) to be converted dtlat,dtlon [ditto] talt Tokyo Datum altitude (in m) to be converted dtalt [ditto] wlat, wlon Calculated results. WGS Lat. and Lon. (in minutes) are given. dwlat,dwlon [ditto] walt Calculated result. WGS altitude (in m) is given. dwalt [ditto] -------------------------------- DRAW Coastline / Lakes and Rivers / Prefecture boundary ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ NAME wshore - Draw coast, lake and rivers with small scale rshore - Define the range to draw coastlines etc. by 'pshore' pshore - Draw coast, lake, rivers, prefecture boundary with larger scale SYNOPSIS call rshore(m, is, in, kw, ke) call rshore(is, in, kw, ke) m Specify whether old Tokyo Datum (m¡á0) or WGS (m¡á1) to use as a reference geodetic system. Must be m=0 or m=1. Otherwise, the calling is treated as the second form with 4 arguments, where old Tokyo Datum (m¡á0) is assumed. is, in Specify latitude range to draw. Must be 20 <= is < in <= 50 . kw, ke Specify longitude range to draw. Must be 120 <= kw < ke <= 155 . external conv call wshore(xorg, yorg, mpen, scl, conv) call pshore(xorg, yorg, mpen, scl, conv) xorg, yorg Coordinates (in cm) of the origin in the drawing. mpen Specify pen numbers to draw coast / lakes and rivers / prefecture boundary in the form of 3 digits number: pen numbers for prefecture boundary / lakes and rivers / coast on the digits of 100 / 10 / 1 , respevtively. As a special case, if mpen=99, the coast of world is drawn by current pen. scl Reduction rate. Coordinate values of scl given by 'conv' is represented by 1 cm of length on the drawing. (If the coordinate values are in km, the drawing scale is 1 : (scl*100,000) .) conv External procedure name to define projection, i.e., the name of subroutine with 4 arguments as follows: subroutine conv(alat, alon, xe, yn) where Lat. (alat) and Lon. (alon) is converted into drawing coordinates toward right (xe) and upward (yn) . SAMPLE external cviken call cvinit(100, 35.*60., 135.*60.) call psopn("wshore.ps","A4P") call plots(2., 0.) call wshore(8., 8., 12, 250., cviken) do 50 i=20,50,5 call cviken(float(i*60), 7200., xe, yn) call plot(8.+xe/250., 8.+yn/250., 3) do 50 k=121,155 call cviken(float(i*60), float(k*60), xe, yn) call plot(8.+xe/250., 8.+yn/250., 2) 50 continue do 60 k=120,155,5 call cviken(1200., float(k*60), xe, yn) call plot(8.+xe/250., 8.+yn/250., 3) do 60 i=21,50 call cviken(float(i*60), float(k*60), xe, yn) call plot(8.+xe/250., 8.+yn/250., 2) 60 continue call plot(0., 16., -3) call wrect(0., 0., 16., 10.) call scisor(0., 0., 16., 10.) call cvinit(253) call cviken(float(35*60), float(135*60), x, y) call rshore(1, 34, 36, 133, 136) call pshore(8.-x/12., 8.-y/12., 212, 12., cviken) call plot(-1., 0., 3) do 70 k=-90,90 call cviken(float(35*60), float(135*60+k), xe, yn) call plot(8.+(xe-x)/12., 8.+(yn-y)/12., 2) 70 continue do 80 k=134,136 call plot(0., -1., 3) do 80 i=-60,30 call cviken(float(35*60+i), float(k*60), xe, yn) call plot(8.+(xe-x)/12., 8.+(yn-y)/12., 2) 80 continue call plote call pscls stop end :/home/SHARE/lib/man/wshore.ps IGRF Calculation ~~~~~~~~~~~~~~~~ NAME igrfs / igrfs2 - Select IGRF model and specify time epoch pgrfs / pgrfs2 - Select PGRF model and specify time epoch dgrfs / dgrfs2 - Select DGRF model and specify time epoch igrfc - Specify position and get total force IGRF value igrfm - Get full components of IGRF vector SYNOPSIS call igrfs(year) call pgrfs(year) call dgrfs(year) call igrfs2(year, mdl) call pgrfs2(year, mdl) call dgrfs2(year, mdl) year Time epoch (in year) of calculation mdl Base year of each model enforced call igrfc(fi, fk, h, f) fi, fk Lat. and Lon. (in degrees) of the point of calculation h Altitude (in m) of the point of calculation f Calculated results. Total force (in nT) value of the model is given. call igrfm(fm) fm Array in which vector values (6 elements) of the model are given. 3 orthogonal components North(X) / East(Y) / Downward (Z) (in nT), and Horizontal(H) component (in nT), Inclination (in degrees) and Declination (in degrees) are put into the array in sequence. USAGE First specify the ModelType (IGRF/DGRF/PGRF) with the Year of Calc. by `call igrfs(year)`, `call igrfs2(year, mdl)`, or another. Then, `call igrfc(FI, FK, H, F)` gives TotalForce (F) of that model at the point of Lat.=FI, Long.=FK, Alt.=H If other components are desired, `call igrfm(FM)` . Here FM is an array with 6 elements, which correspond to North(X), East(Y), Downward(Z), Horizontal(H) components, Inclination(I) and Declination(D). -------------------------------- Synthetic Magnetic Anomaly of Trivial model ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ NAME magafd - Define magnetization and ambient magnetic field direction mpoint - Specify point source and its parameter mvline - Specify vertical line source and its parameter mhrect - Specify horizontal rectangle surface source and its parameter mprism - Specify rectangular block source and its parameter calma - Calculate magnetic anomaly at specified observation point (Note) Assume righthanded Cartesian coordinates X pointing North, Y pointing East, and Z pointing Down. All units of lengths can be arbitrary as far as they are common to all lengths. SYNOPSIS call magafd(am, the, phi, dip, dec) am Magnetization intensity (in A/m) the, phi Inclination and declination (in degrees) of magnetization dip, dec Inclination and declination (in degrees) of ambient field (Earth's magnetic field) call mpoint(xc, yc, zc, vol) call mvline(xc, yc, zs, w, area) call mhrect(xs, xt, ys, yt, zc, h) call mprism(xs, xt, ys, yt, zs, w) xc, yc Horizontal coordinates of source center zc Downward coordinate of source center xs, xt X-axis range (lower and upper limits) of source ys, yt Y-axis range (lower and upper limits) of source zs, w Top deapth and thickness in depth range If w=0., bottom unlimitted model is assumed. vol Volume represented by a point (in the unit of cube of common length unit) area Area of cross-section represented by a line (in the unit of square of common length unit) h Thickness of the body represented by a surface (in the unit common with lengths) f = calma(xp, yp, zp) xp, yp Horizontal coordinates of the point of calculation zp Depth coordinate of the point of calculation f Calculated result. Total force (in nT) value of magnetic anomaly is given. -------------------------------- Regression Analysis ~~~~~~~~~~~~~~~~~~~ NAME sm1opn / sm2opn / sm3opn - Initialize regression analysis sm1ex / sm2ex / sm3ex - Indicate each data to apply regression analysis sm1cls / sm2cls / sm3cls - Calculate regression coefficients sm1rv / sm2rv / sm3rv - Calculate the value(s) of regression formula sm1opn / sm1ex / sm1cls / sm1rv : for N dimension linear regression (N <= 8) sm2opn / sm2ex / sm2cls / sm2rv : for 2-dimensional N-th degree regression (N <= 8) sm3opn / sm3ex / sm3cls / sm3rv : for 3-dimensional N-th degree regression (N <= 5) for 5 dimension (t1-t5) linear regression: c0 + c1*t1 + c2*t2 + c3*t3 + c4*t4 + c5*t5 for 2 dimension (u,v) 4th degree regression: c0 + c10*u + c11*v + c20*u*u + c21*u*v + c22*v*v + c30*u*u*u + c31*u*u*v + c32*u*v*v + c33*v*v*v + c40*u*u*u*u + c41*u*u*u*v + c42*u*u*v*v + c43*u*v*v*v + c44*v*v*v*v for 3 dimension (x,y,z) cubic regression: c0 + c11*x + c12*y + c13*z + c211*x*x + c212*x*y + c213*y*z + c222*y*y + c223*y*z + c233*z*z + c3111*x*x*x + c3112*x*x*y + c3113*x*x*z + c3122*x*y*y + c3123*x*y*z + c3133*x*z*z + c3222*y*y*y + c3223*y*y*z + c3233*y*z*z + c3333*z*z*z SYNOPSIS call sm1opn(nvar, nsrc1) call sm2opn(ndg2, nsrc2) call sm3opn(ndg3, nsrc3) nvar Number of dimension (1 <= nvar <= 8) ndg2 Number of degree for 2-dimensional regression (1 <= ndg2 <= 8) ndg3 Number of degree for 3-dimensional regression (1 <= ndg3 <= 5) nsrc [nsrc1, nsrc2, nsrc3] Number of variables to apply regression (1 <= nsrc <= 5) call sm1ex(var, src1) call sm2ex(u, v, src2) call sm3ex(x, y, z, src3) var Array containing independent parameters (nvar elements) u, v 2 independent parameters of 2-dimensional regression x, y, z 3 independent parameters of 3-dimensional regression src [src1, src2, src3] Array containing variables to apply regression (nsrc elements) call sm1cls(coef, m, n) call sm2cls(coef, m, n) call sm3cls(coef, m, n) coef Array to return regression coefficients m, n Adjustable size of coef [dimension coef(m,n)] Usually, m and n are set to the number of regression coefficients and the number of variables (nsrc). If m and/or n is less than those numbers, the excess part of coefficients are not returned. So, m and/or n may be zero or minus. Even if m and/or n is greater than those numbers, only actual numbers of coefficients are returned, and further location of the array is not changed. Number of regression coefficients are: 2, 3, 4, 5, 6, 7, 8, 9 for N dimension linear regression (N=1,2,3,4,5,6,7,8) , 3, 6,10,15,21,28,36,45 for 2-dimensional N-th degree regression (N=1,2,3,4,5,6,7,8) , 4,10,20,35,56 3-dimensional N-th degree regression (N=1,2,3,4,5) . call sm1rv(var, sml) call sm2rv(u, v, sml) call sm3rv(x, y, z, sml) var Array defining the independent parameters to calculate regression formula value(s). (nvar elements) u, v 2 independent parameters to calculate 2-dimensional regression formula value(s). x, y, z 3 independent parameters to calculate 3-dimensional regression formula value(s). sml Array to return regression formula value(s) (nsrc elements) -------------------------------- Random Number Generation ~~~~~~~~~~~~~~~~~~~~~~~~ NAME rand1 - Generate uniform random number within the range of [-1., +1.] randg - Generate Gaussian random number with zero mean and dispersion 1. SYNOPSIS rv1 = rand1() rvg = randg() SOURCE LIST #include #include double rand1(void) { return(drand48()*2.-1.); } double randg(void) /* after Box-Muller method */ { static int k=0; static double u, v, r, f; while (k == 0) { u = rand1(); v = rand1(); r = u*u+v*v; if ((r<1.) && (r>0.)) { f = sqrt(-2.*log(r)/r); k = 1; return(u*f); } } k = 0; return(v*f); } float rand1_(void) { return((float)rand1()); } float randg_(void) { return((float)randg()); } -------------------------------- Assistance to Show the progress of loop operation, ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ to Set-up Working Directory Path, ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ and to Read-in Process Parameters ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ NAME premsg / prompt - Output message to screen with/without NewLine dpcent - Display the percentage of the progress of loop operation. opnpin / clspin - Provide a chance to select predefined data file as a data source for reading process parameters lwkdir - Display HOME (~) and CURRENT (.) directory paths, prompt and wait for specifying Working Directory Path, and then dispay absolute path name of that directory and filename list in it. The absolute pathname of the directory is returned as a string (with trailing '/'), which can be used to specify absolute filename simply by appending short filename. Here the "Working Directory" does not mean the change of Current Working Directory environment, but offers the easier way of specifying absolute path to target file. parmin / gparma / gparmi / gparmf / gparmd / gparmif / gparmid / gparmi2 / gparmf2 / gparmd2 - Output prompt message and read process parameter(s) SYNOPSIS call premsg(mesg) call prompt(mesg) mesg Message string to output Output message to screen ("stderr" output) with NewLine ("premsg") or without NewLine ("prompt") call dpcent(m, n) n Total number of loops. m Number of loops executed already. At the first call, give m=0, n=0 to initialize. Then call it each time m=1,2,...,n (n>0). call opnpin() Give a chance to select predefined data file for parameter input. If NULL string is specified, process parameters are read from Keyboard as usual. This mechanism is effective for "lwkdir" function and "parmin/gparm?i/gparm??" routines until "clspin" is called. call clspin() Finalize the "opnpin" mechanism. nch = lwkdir(len, dnm) len Maximum length of string variable 'dnm'. dnm String variable to return the absolute pathname of Working Directory. (must be of size 'len' or more) nch Number of characters returned in 'dnm' The string returned in 'dnm' is in the form of C language string having a trailing NULL (char(0)) character. If the length of path string of specified working directory is equal to or higher than 'len', the process is terminated abnomally with an error message output. If the length of path string of HOME (~), CURRENT (.), or Working Directory is higher than 120, the process is terminated abnomally with an error message output. call parmin(str, lbuf, buf) str Prompt message string lbuf Maximum size (in bytes) to input to buf buf Raw input string Output prompt message and input raw string under "opnpin" control. call gparma(str, lnam, nam) call gparmi(str, iv) call gparmf(str, fv) call gparmd(str, dv) call gparmif(str, iv, fv) call gparmid(str, iv, dv) call gparmi2(str, iv, jv) call gparmf2(str, fv, gv) call gparmd2(str, dv, ev) str Prompt message string lnam Maximum size (in bytes) to input to nam nam String parameter without in-between white space is given. iv,jv Integer value parameter is given. fv,gv Floating value parameter is given. dv,ev Double floating value parameter is given. Output prompt message and read parameter(s) under "opnpin" control. EXAMPLE character wdr*80, fnam*120 c call premsg('--- start ---') call opnpin() ldr = lwkdir(80,wdr) + 1 fnam = wdr call gparma('Enter filename ==> ', 121-ldr, fnam(ldr:120)) open(10,file=fnam,status='old') : : call clspin() call prompt('Loop operation - ') call dpcent(0, 0) do 10 i=1,imax do 10 j=1,jmax : : call dpcent((i-1)*jmax+j, imax*jmax) 10 continue : :