program spheresolver
!c Combined program to find Weber solution in 3d space and project that
!c solution to the surface of a unit sphere, and to output objective function
!c data over the projection point area. Maximum number of demand points is 100.
real x,y,z,x1,x2,x3,y1,y2,y3,c,c1,c2,c3,v1,v2,v3,pi,deltax,deltay,deltaz,xstart,ystart,zstart,error
real latfinal,longfinal,tempsum,maximumdistance,centertheta,centerphi,radius,solutionvalue,lowest
real istep,jstep,highest,gammasum,fvalue
integer howmany,iteration,i,j,k,maxpoint1,maxpoint2,flag,matches,lowestxposition,lowestyposition
integer highestxposition,highestyposition
parameter (pi=3.141592654)
dimension phi(100),theta(100),weight(100),xdata(100),ydata(100),zdata(100),distance(100,100)
dimension sumdst(360,360),gamma(100),greek(100),temp(100)
2 format (f12.8,',',f12.8,',',f12.8)
3 format (f5.2,',',f5.2,',',f12.8)
print*
print*
print*
print*, 'WELCOME TO THE 1-MEDIAN SPHERICAL SURFACE PROBLEM SOLVER '
print*, '2001, Impact Consulting. Unlimited distribution.'
print*, 'This program solves a 3-dimensional Weber problem for user-defined'
print*, 'demand points on the surface of a unit sphere, projects the solution'
print*, 'to the surface of the sphere, and returns diagnostic data regarding'
print*, 'solution optimality. Lat/long are based on the program convention.'
print*, 'From 0 Lat 0 Long (on the prime meridian at the equator), proceed'
print*, 'Lat positive headed north and negative headed south, and proceed Long'
print*, 'positive headed east and negative headed west.'
print*
!c user defines how many demand points will be analyzed
print*, 'Enter how many demand points will be analyzed.'
read*, howmany
print*
!c user enters demand point data for (howmany) arrays.
!c theta and phi data are converted to radians and stored in arrays
do 100 i=1,howmany
print*, 'DEMAND POINT',i
print*, 'Enter longitude (theta) for point number',i,'.'
read*, theta(i)
!c convert longitude to radians
theta(i)=theta(i)*pi/180
print*, 'Enter latitude for point number',i,'.'
read*, temp(i)
!c convert latitude to radians then to colatitude
temp(i)=temp(i)*pi/180
phi(i)=(pi/2)-temp(i)
print*, 'Enter positive weight for point number',i,'.'
read*, weight(i)
100 continue
!c program now creates array of x,y,z positions for the data points
!c for the unit sphere (r=1)
do 110 i=1,howmany
xdata(i)=cos(theta(i))*sin(phi(i))
ydata(i)=sin(theta(i))*sin(phi(i))
zdata(i)=cos(phi(i))
110 continue
do 111 i=1,howmany
!c system check - print data to screen for all demand points in deg and xyz coords
print*, 'point',i,'long/colat/w',(theta(i)*180/pi),(phi(i)*180/pi),weight(i)
print*, 'x,y,z',xdata(i),ydata(i),zdata(i)
111 continue
!c program solves the 3d Weber problem based on the x,y,z data for
!c each demand point using Weiszfeld's algorithm and a small error
xstart=0
ystart=0
zstart=0
xnew=0
ynew=0
znew=0
deltax=1
deltay=1
deltaz=1
iteration=0
error=.00000001
120 if ((deltax.gt.(0.00000001)).or.(deltay.gt.(0.00000001)).or.(deltaz.gt.(0.00000001))) then
iteration=iteration+1
do 130 i=1,howmany
gamma(i)=(weight(i)/(((xstart-xdata(i))**2+(ystart-ydata(i))**2+(zstart-zdata(i))**2+error)**(0.5)))
130 continue
gammasum=0
do 135 i=1,howmany
gammasum=gammasum+gamma(i)
135 continue
do 140 i=1,howmany
greek(i)=(gamma(i)/gammasum)
140 continue
xnew=0
ynew=0
znew=0
do 150 i=1,howmany
xnew=xnew+(greek(i)*xdata(i))
ynew=ynew+(greek(i)*ydata(i))
znew=znew+(greek(i)*zdata(i))
150 continue
!c compute the function value (minisum) for the 3d L2 problem and report results
! fvalue=0
! do 155 i=1,howmany
! fvalue=fvalue+ weight(i)*sqrt((xdata(i)-xnew)**2+(ydata(i)-ynew)**2+(zdata(i)-znew)**2)
!155 continue
! print*, 'iteration',iteration,'3d minisum value=,',fvalue
! print*, 'at x,y,z=',xnew,ynew,znew
! read*
deltax=abs(xstart-xnew)
deltay=abs(ystart-ynew)
deltaz=abs(zstart-znew)
xstart=xnew
ystart=ynew
zstart=znew
go to 120
endif
!c print number of iterations required to arrive at solution
print*
print*, 'Problem converged in',iteration,'iterations.'
!c print x,y,z solution
print*, '3d Weber solution x,y,z coordinates are'
print 2, xnew,ynew,znew
print*
read*
!c print latitude,longitude solution
print*, 'The projected solution point on the sphere is located at'
longfinal=projectlongitude(xnew,ynew)
latfinal=projectlatitude(xnew,ynew,znew)
if (latfinal.lt.0) then
latfinal=latfinal+pi
endif
print*, longfinal*180/pi,'longitude,'
print*, latfinal*180/pi,'colatitude, and'
print*, ((pi/2)-latfinal)*180/pi,'latitude.'
do 50 i=1,howmany
phi(i)=temp(i)
50 continue
maximumdistance=0
!c program determines the smallest spherical circle containing the demand points
!c iterates through each point combination to find the maximum distance
!c searches and retains a greatest i,j combination
do 180 i=1,howmany
do 170 j=1,howmany
distance(i,j)=arcdistance(theta(i),phi(i),theta(j),phi(j))
! print*, theta(i),phi(i),theta(j),phi(j)
! print*, i,j,distance(i,j)*180/pi
if (distance(i,j).gt.maximumdistance) then
maximumdistance=distance(i,j)
maxpoint1=i
maxpoint2=j
endif
170 continue
180 continue
!c determines spherical circle minimum radius
radius=maximumdistance/2
if (radius.gt.(pi/4)) then
print*, 'The spherical radius is larger than the pi/4 restriction ensuring'
print*, 'optimality. The absolute minimum radius based on farthest points is'
print*, radius/pi*180,'degrees, although the spherical covering circle may be larger.'
else
print*, 'The spherical radius may (but is not guaranteed to) be within the'
print*, 'pi/4 restriction ensuring optimality. The minimum radius based on'
print*, 'farthest points is',radius/pi*180,'degrees, although the covering circle'
print*, 'may be larger.'
endif
print*, 'Two points defining a possible absolute minimum covering circle are'
print*, 'demand points',maxpoint1,'and',maxpoint2,'. Other points with the same'
print*, 'distance as these points may exist and these points are not necessarily'
print*, 'on the spherical circle edge.'
!c compute the objective function value at the projected point
solutionvalue=0
do 3000 k=1,howmany
solutionvalue=solutionvalue+weight(k)*arcdistance(theta(k),temp(k),longfinal,((pi/2)-latfinal))
3000 continue
print*, 'The projected point objective function value is',solutionvalue,'.'
!c the program now takes divisions across the sphere to compute
!c objective function values for plotting (degree by degree plot) which will be
!c heavily weighted towards the poles (more data computed there)
!c compute each demand point's weighted distance to the gridpoint intersection
open (13,file='spherex.loc')
write(13,*) 'longitude, latitude, objfuncvalue'
print*, 'Writing to file - please be patient.'
!c compute/write for positive (E) long, positive (N) lat
do 2000 i=1,90
do 2005 j=1,44
sumdst(i,j)=0
2005 continue
2000 continue
do 2010 i=1,90
do 2020 j=1,44
do 2030 k=1,howmany
gridpointdistance=arcdistance(theta(k),temp(k),((i-1)*pi/90),(j*pi/90))
gridpointdistance=abs(gridpointdistance*weight(k))
sumdst(i,j)=sumdst(i,j)+gridpointdistance
2030 continue
write(13,*) (i-1)*2,',',j*2,',',sumdst(i,j)
2020 continue
2010 continue
!c compute/write for negative (W) long, positive (N) lat
do 2040 i=1,90
do 2045 j=1,44
sumdst(i,j)=0
2045 continue
2040 continue
do 2050 i=1,90
do 2060 j=1,44
do 2070 k=1,howmany
gridpointdistance=arcdistance(theta(k),temp(k),(i*pi/90*(-1)),(j*pi/90))
gridpointdistance=abs(gridpointdistance*weight(k))
sumdst(i,j)=sumdst(i,j)+gridpointdistance
2070 continue
write(13,*) (i*(-1))*2,',',j*2,',',sumdst(i,j)
2060 continue
2050 continue
!c compute/write for positive (E) long, negative (S) lat
do 2071 i=1,90
do 2075 j=1,44
sumdst(i,j)=0
2075 continue
2071 continue
do 2080 i=1,90
do 2090 j=1,44
do 2100 k=1,howmany
gridpointdistance=arcdistance(theta(k),temp(k),((i-1)*pi/90),(j*pi/90*(-1)))
gridpointdistance=abs(gridpointdistance*weight(k))
sumdst(i,j)=sumdst(i,j)+gridpointdistance
2100 continue
write(13,*) (i-1)*2,',',(j*(-1))*2,',',sumdst(i,j)
2090 continue
2080 continue
!c compute/write for negative (W) long, negative (S) lat
do 2110 i=1,90
do 2115 j=1,44
sumdst(i,j)=0
2115 continue
2110 continue
do 2120 i=1,90
do 2130 j=1,44
do 2140 k=1,howmany
gridpointdistance=arcdistance(theta(k),temp(k),(i*pi/90*(-1)),(j*pi/90*(-1)))
gridpointdistance=abs(gridpointdistance*weight(k))
sumdst(i,j)=sumdst(i,j)+gridpointdistance
2140 continue
write(13,*) (i*(-1))*2,',',(j*(-1))*2,',',sumdst(i,j)
2130 continue
2120 continue
!c compute/write for poles and projected point
do 2150 k=1,howmany
gridpointdistance=arcdistance(theta(k),temp(k),0,90)
gridpointdistance=abs(gridpointdistance*weight(k))
sumdst(i,j)=sumdst(i,j)+gridpointdistance
2150 continue
write(13,*) 0,',',90,',',sumdst(i,j)
do 2160 k=1,howmany
gridpointdistance=arcdistance(theta(k),temp(k),0,-90)
gridpointdistance=abs(gridpointdistance*weight(k))
sumdst(i,j)=sumdst(i,j)+gridpointdistance
2160 continue
write(13,*) 0,',',-90,',',sumdst(i,j)
write(13,*) longfinal*180/pi,',',((pi/2)-latfinal)*180/pi,',',solutionvalue
close (13)
print*, 'Data has been written to file for external processing.'
print*, 'Press enter to exit the program.'
read*
!c end program
end
function projectlatitude (x,y,z) result (m)
real m
m=atan((sqrt(x**2+y**2)/z))
continue
end function projectlatitude
function projectlongitude (x,y) result (n)
real n
n=atan(y/x)
continue
end function projectlongitude
function arcdistance (theta1,phi1,theta2,phi2) result (p)
real p
p=acos((cos(phi1)*cos(phi2)*cos(theta1-theta2))+(sin(phi1)*sin(phi2)))
continue
end function arcdistance