Reputation: 1097
Project description:
I have the coordinates of McDonalds in cities and the coordinates of city centers. They are divided to text files by city. For example,
=== atlanta_mc.txt === === houston_mc.txt === ... (for 200 cities)
[latitude] [longitude] [latitude] [longitude] <-- The headers are not written in the files.
33.5431140 -84.5766910 29.8044570 -95.3990780
34.0489350 -84.0928960 30.0051170 -95.2834550
33.7853660 -84.1018980 29.7371140 -95.5352550
34.0903210 -84.2040180 29.7280160 -95.4188230
34.1606520 -84.1781010 29.7851700 -95.3609770
.... ....
My goal is to calculate the distance between the city center and each McDonald for each city. The city center data file has [cityname],[latitude],[longitude] fiels. It looks like:
== city_centers.txt ==
[city name] [latitude] [longitude] <-- The headers are not written in the file.
atlanta 33.79526 -84.326528
la 34.057453 -118.413842
ny 40.759347 -73.980202
houston 29.733215 -95.430824
....
(for 200 cities in the US.)
Fortran code:
For one city (say, Atlanta), I can do the job as follows:
program mcdonalds
implicit none
interface
function distkm(deglat1,deglon1,deglat2,deglon2)
real :: distkm
real, intent(in) :: deglat1,deglon1,deglat2,deglon2
end function distkm
end interface
integer, parameter :: maxnr = 200
integer :: nr, i, j, ios
character(len=1) :: junkfornr
real, dimension(:), allocatable :: lat, lon
! Obtain the total number of observations
nr=0
open(10, file='atlanta_mc.txt', status='old')
do i=1,maxnr
read(10,*,iostat=ios) junkfornr
if (ios/=0) exit
if (i == maxnr) then
stop
endif
nr = nr + 1
end do
allocate(lat(nr))
allocate(lon(nr))
rewind(10)
do i=1,nr
read(10,*) lat(i), lon(i)
end do
open(20,file='RESULT_atlanta_mc.txt')
do i=1,nr
write(20,100) "atlanta",lat(i),lon(i),distkm(lat(i),lon(i),33.79526,-84.326528)
100 format(A,3F12.6)
end do
end program mcdonalds
function distkm(deglat1,deglon1,deglat2,deglon2)
implicit none
real, intent(in) :: deglat1,deglon1,deglat2,deglon2
real :: pi,dlat,dlon,lat1,lat2,a,distkm
pi = 4*atan(1.0)
dlat = (deglat2-deglat1)*pi/180
dlon = (deglon2-deglon1)*pi/180
lat1 = deglat1*pi/180
lat2 = deglat2*pi/180
a = (sin(dlat/2))**2 + cos(lat1)*cos(lat2)*(sin(dlon/2))**2
distkm = 6372.8*2*atan2(sqrt(a),sqrt(1-a))
end function distkm
Problem:
I have to do this for the all cities. So, I have to change (1) the McDonalds data file, (2) the city center coordinate; and (3) the result file for each city accordingly. How do I have to write a program to iterate this process for the whole 200 cities? I use Intel Fortran (ifort) to compile it.
EDIT (This code worked -- Friday, April 11, 2014):
! calculate distances from the city center to each McDonald
! Try this for all the cities
program mcdonalds
implicit none
interface
function distkm(deglat1,deglon1,deglat2,deglon2)
real :: distkm
real, intent(in) :: deglat1,deglon1,deglat2,deglon2
end function distkm
end interface
integer, parameter :: maxnr = 200
integer :: nr, i, j, ios
character(len=1) :: junkfornr
character(len=25) :: fn
real, dimension(:), allocatable :: lat, lon
real, dimension(5) :: clat, clon
character(len=15), dimension(5) :: filename
character(len=6), dimension(5) :: code
character(len=7), dimension(5) :: cityname
! I have a text file that contains the list of data file names, city names and the (lat/lon)s of city centers.
! The name of that file: datafilelistcitycenter.txt
! === The content of the file looks like: ===
! atlanta_mc.txt 11111 atlanta 33.79526 -84.326528
! houston_mc.txt 22222 houston 29.733215 -95.430824
! la_mc.txt 33333 la 34.057453 -118.413842
! ny_mc.txt 44444 ny 40.759347 -73.980202
! philly_mc.txt 55555 philly 29.733215 -95.430824
! Let's assume we know the number of the cities.
! In this case, it is 5.
open(10, file='datafilelistcitycenter.txt', status='old')
rewind(10)
do i=1,5
read(10,*) filename(i), cityname(i), clat(i), clon(i)
end do
close(10)
! Obtain the total number of observations for each city
do j=1,5
open(j, file=filename(j), status='old')
nr=0
do i=1,maxnr
read(j,*,iostat=ios) junkfornr
if (ios/=0) exit
if (i == maxnr) then
write(*,*) "Error: Maximum number of records exceeded..."
write(*,*) "Exiting the program..."
stop
endif
nr = nr + 1
end do
allocate(lat(nr),lon(nr))
rewind(j)
do i=1,nr
read(j,*) lat(i), lon(i)
end do
close(j)
write(fn, '("RESULT_",A15)') filename(j)
open(j,file=fn)
do i=1,nr
write(j,100) code(j),cityname(j),lat(i),lon(i),distkm(lat(i),lon(i),clat(j),clon(j))
100 format(A6,A7,3F12.6)
end do
close(j)
deallocate(lat,lon)
end do
end program mcdonalds
function distkm(deglat1,deglon1,deglat2,deglon2)
implicit none
real, intent(in) :: deglat1,deglon1,deglat2,deglon2
real :: pi,dlat,dlon,lat1,lat2,a,distkm
pi = 4*atan(1.0)
dlat = (deglat2-deglat1)*pi/180
dlon = (deglon2-deglon1)*pi/180
lat1 = deglat1*pi/180
lat2 = deglat2*pi/180
a = (sin(dlat/2))**2 + cos(lat1)*cos(lat2)*(sin(dlon/2))**2
distkm = 6372.8*2*atan2(sqrt(a),sqrt(1-a))
end function distkm
No luck. I end up with very weird messages: severe <66>: output statement overflows record, unit -5, file Internal Fomatted Write.
EDIT: I think I got it right. Thanks.
Upvotes: 0
Views: 146
Reputation: 29401
Instead of opening one specific city file such as atlanta_mc.tx, create a file that a list of the filenames of the city. Open that file, add a new loop to your program and read a name of a city file, open and process it as you do now, then close the file. Next iteration of new loop will process the next city file, etc. This is assuming that the cities are in separate files. If they are in one big file, then create arrays to hold information for all of the cities, read all the information into the array, then have a new loop to iterate over the cities.
Also, if you place your function distkm
into a module you won't have to bother with having an interface
for it. Using modules and easier and less error prone. An example: distkm
Upvotes: 2