program greatcircle
!c Program takes latitude, longitude, and w inputs for problems and atyts to
!c find solutions to Weiszfeld's algorithm using the great circle arc distance metric.
!c This is not expected to be highly successful, but the results need to be compared to
!c those from the developed geometric methods.
real pi,deltax,deltay,xstart,ystart,ysum,xnew,ynew,error
integer howmany,i,j,k
parameter (pi=3.141592654)
dimension x(100),y(100),w(100),a(100),b(100)
2 format (f12.8,',',f12.8,',',f12.8)
3 format (f5.2,',',f5.2,',',f12.8)
xstart=-45*pi/180
ystart=-45*pi/180
deltax=1
deltay=1
error=.000001
print*
print*
print*
print*, 'WELCOME TO THE WEISZFELD GREAT CIRCLE ARC DISTANCE TOOL'
print*, '2001, Impact Consulting. Unlimited distribution.'
print*, 'This program atyts to use the Weiszfeld algorithm and a great circle'
print*, 'arc norm to find a Weber problem solution on the surface of a unit'
print*, 'sphere. It returns a point in lat/long for comparison with solutions'
print*, 'from other methods. 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 50 i=1,howmany
print*, 'DEMAND POINT',i
print*, 'Enter longitude (theta) for point number',i,'.'
read*, x(i)
!c convert longitude to radians
x(i)=x(i)*pi/180
print*, 'Enter latitude for point number',i,'.'
read*, y(i)
!c convert latitude to radians then to colatitude
y(i)=y(i)*pi/180
print*, 'Enter positive weight for point number',i,'.'
read*, w(i)
50 continue
100 if (deltax.gt.(0.0001).or.deltay.gt.(0.0001)) then
do 110 i=1,howmany
a(i)=w(i)/(arcdistance(xstart,ystart,x(i),y(i))+error)
110 continue
ysum=0
do 120 i=1,howmany
ysum=ysum+a(i)
120 continue
do 130 i=1,howmany
b(i)=a(i)/ysum
130 continue
xnew=0
ynew=0
do 200 i=1,howmany
xnew=xnew+(b(i)*x(i))
if (xnew.gt.180) then
xnew=(360-xnew)*(-1)
endif
ynew=ynew+(b(i)*y(i))
if (ynew.gt.90) then
ynew=ynew-90
xnew=(360-xnew)*(-1)
endif
if (ynew.lt.(-90)) then
ynew=-90-(ynew+90)
xnew=(360-xnew)*(-1)
endif
200 continue
print*, 'iteration x,y resultant is',xnew*180/pi,ynew*180/pi
deltax=abs(xnew-xstart)
deltay=abs(ynew-ystart)
xstart=xnew
ystart=ynew
go to 100
endif
print*, 'Press enter to exit the program.'
read*
!c end program
end
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