subroutine mpstrt(stcprm, conang, p_lat, p_long, r_lat, r_long) parameter (pi=3.14159265358979323846,radpdg=pi/180.) parameter (dgprad=180./pi) parameter (REARTH=6371.2) dimension stcprm(15) Ctypedef struct { C 1: gamma, C 2-10: rotate(3,3), C 11-12: x0,y0, C 13-14: crotate,srotate, C 15: grdszq C C* C* General Purpose routine to Set Map Parameters. Called by the C* special purpose routines for oblique stereographic, oblique and C* transverse Mercator, oblique Lambert Conformal, etc. Projections C* Inputs: p_lat,p_long - Latitude and Longitude of the Pole C* Point for the Projection C* r_lat,r_long - Latitude and Longitude of a C* reference point - 180degrees around the Pole Point C* from the cut C* conang - angle between the Projection Pole axis and C* the generators (sides) of the cone on which the Earth C* is projected. + or - 90 degrees indicates a plane, C* hence Stereographic; 0 degrees indicates a cylinder C* and hence Mercator. C* Outputs: stcprm - map parameters C*/ dimension temp(3) ifind(l,k) = -2 + 3*l + k call ll_geo(p_lat,p_long,temp) do k=1,3 stcprm(ifind(3,k)) = temp(k) enddo C for (k=0;k<3;k++) stcprm->rotate[2][k] = temp.v[k]; call ll_geo(r_lat,r_long,temp) do k=1,3 stcprm(ifind(1,k)) = temp(k) enddo C for (k=0;k<3;k++) stcprm->rotate[0][k] = temp.v[k]; tnorm = x_prod(stcprm(ifind(3,1)),stcprm(ifind(1,1)), C stcprm(ifind(2,1))) C norm = C x_product (stcprm->rotate[2],stcprm->rotate[0],stcprm->rotate[1]); do k=1,3 stcprm(ifind(2,k)) = stcprm(ifind(2,k)) /tnorm enddo C for (k=0;k<3;k++) stcprm->rotate[1][k] /= norm tnorm = x_prod(stcprm(ifind(2,1)),stcprm(ifind(3,1)), C stcprm(ifind(1,1))) C x_product (stcprm->rotate[1],stcprm->rotate[2],stcprm->rotate[0]); C stcprm->x0 = stcprm->y0 = stcprm->srotate = 0; stcprm(11)=0. stcprm(12)=0. stcprm(14)=0. stcprm(13)=1. stcprm(15) = REARTH stcprm(1) = sin(RADPDG * conang ) C stcprm->crotate = 1.; C stcprm->gridszeq = REARTH; C stcprm->gamma = sin(RADPDEG * cone_ang); C*geographic triple : i = equator @ std merid, C j = equator @ 90E, C k = North Pole */ C*map base triple : i' = M-prime meridian @ M-equator, C j' = M-East at location i', C k' = M-Pole */ C return end