http://www.aiuxian.com/article/p-2767848.html
#define toDeg(X) (X*180.0/M_PI)
/**
* @method 根据两点经纬度，计算与真北方夹角
*
* @param longitude1
* @param latitude1
* @param longitude2
* @param latitude2 // 目标点
*/
- (void)getAngle:(double)longitude1 latitude1:(double)latitude1 longitude2:(double)longitude2 latitude2:(double)latitude2 {
double cos_c = cos(90 - latitude2)*cos(90 - latitude1) + sin(90 - latitude2)*sin(90 - latitude1)*cos(longitude2-longitude1);
double sin_c = sqrt(1 - pow(cos_c, 2));
double z = asin(sin(90 - latitude2)*sin(longitude2 - longitude1)/sin_c);
z = toDeg(z);
// A（起始点）为原点B目标点
if(longitude1 < longitude2 && latitude1 < latitude2) { // 第一象限
} else if (longitude1 < longitude2 && latitude1 > latitude2){ // 第二象限
z += 360;
} else { // 三四象限
z = 180 - z;
}
cityHeading = z;
NSLog(@"城市夹角:%f",z);
}
此处设定求B相对于A的方位角，即A为当前位置，B为目标位置
Aj：A点经度
Aw：A点纬度
Bj：B点经度
Bw：B点纬度
北纬为正，南纬为负；东经为正，西经为负
经纬度使用度，DDD.DDDDDD°，非度分或度分秒。
度数未加说明均采用角度制
R：地球平均半径
Azimuth：方位角，以真北为0度起点，由东向南向西顺时针旋转360度
这里需要注意一点，我们一开始的假设便是求B点相对于A点的方位角，因此这里是Bj-Aj，不要写反，否则得不到正确结果。
算到这里，还没有完，得到的结果并不总符合我们对方位角的定义，因此要根据B相对于A的位置在四个象限两个轴上进行讨论，依据不同情况对计算结果进行不同处理。假设A点固定于原点，则:
B点在第一象限，Azimuth=A；
B在第二象限，Azimuth=360+A;
B在第三四象限，Azimuth=180-A。
这里只说了象限的讨论结果，因为轴上的讨论更复杂些，要结合程序运行环境一起考虑，考虑的主要因素是系统的计算精度。譬如，在三面角余弦公式中，当AB点纬度值相同时，对公式的值起决定作用的就是cos(Bj-Aj)这一项，当Bj-Aj的值比较小时，例如0.0001（这在赤道地区对应的长度为11米左右），用一般的计算器计算时值为1，这样，后面的计算便不可能完成。但是，如果用计算机计算则为0.999999999998476913…………。所以，基于以上原因，需要对轴的“范围进项扩充”，要用单片机、手机运算的尤其要注意。
经过一系列计算，最后，就得到了最终结果。
似乎有人注意到了，以上的计算都是把地球看成标准的球体，而事实是地球是个椭圆，其实，地球的偏心率极低，各位可以将此法得到的计算结果与谷歌地球（WGS84坐标系统，我说的不是谷歌地图）上的结果进行对比，偏差是非常小的（我测的几个值，最大偏差0.5度）。
二、距离的求算
其实，“眼尖”的或许已经注意到了，第一步的余弦值结果就可以直接用来求算AB两点间的球面距离，用反余弦函数求得c的度数，再将度数转换为弧度，乘以地球半径就得到了两点间的球面距离。
公式为
这里要注意，L的单位与R的单位一致，单位不同的不要忘记换算
短距离（例如100米，30米）使用这个公式，计算出的结果与谷歌地球给出的距离偏差在0.5%以下，长距离计算时，偏差则可以降至0.01%以下。求算的距离越大，偏差越小，就是这个公式的特点，原因不说自明。
PS:对于一些GPS接收机，其数据格式为NMEA-0183，经纬度数据为DDDMM.MMMM，需要将它转换为度，公式为：
经纬度（度）=DDD+MM.MMMM/60
第二种方法：
1 #import <Foundation/Foundation.h> 2 3 @interface MyLatLng : NSObject { 4 double m_LoDeg,m_LoMin,m_LoSec; 5 double m_LaDeg,m_LaMin,m_LaSec; 6 double m_Longitude,m_Latitude; 7 double m_RadLo,m_RadLa; 8 double Ec; 9 double Ed; 10 } 11 - (id)init:(double)longitude latitude:(double)latitude; 12 13 @property (assign, nonatomic) double m_LoDeg; 14 @property (assign, nonatomic) double m_LoMin; 15 @property (assign, nonatomic) double m_LoSec; 16 @property (assign, nonatomic) double m_LaDeg; 17 @property (assign, nonatomic) double m_LaMin; 18 @property (assign, nonatomic) double m_LaSec; 19 @property (assign, nonatomic) double m_Longitude; 20 @property (assign, nonatomic) double m_Latitude; 21 @property (assign, nonatomic) double m_RadLo; 22 @property (assign, nonatomic) double m_RadLa; 23 @property (assign, nonatomic) double Ec; 24 @property (assign, nonatomic) double Ed; 25 @end
1 #import "MyLatLng.h" 2 #define RC 6378137 3 #define RJ 6356725 4 5 @implementation MyLatLng 6 @synthesize m_LoDeg; 7 @synthesize m_LoMin; 8 @synthesize m_LoSec; 9 @synthesize m_LaDeg; 10 @synthesize m_LaMin; 11 @synthesize m_LaSec; 12 @synthesize m_Longitude; 13 @synthesize m_Latitude; 14 @synthesize m_RadLa; 15 @synthesize m_RadLo; 16 @synthesize Ec; 17 @synthesize Ed; 18 19 - (id)init:(double)longitude latitude:(double)latitude{ 20 self = [super init]; 21 22 if (self) { 23 m_LoDeg=(int)longitude; 24 m_LoMin=(int)((longitude-m_LoDeg)*60); 25 m_LoSec=(longitude-m_LoDeg-m_LoMin/60.)*3600; 26 27 m_LaDeg=(int)latitude; 28 m_LaMin=(int)((latitude-m_LaDeg)*60); 29 m_LaSec=(latitude-m_LaDeg-m_LaMin/60.)*3600; 30 31 m_Longitude=longitude; 32 m_Latitude=latitude; 33 m_RadLo=longitude*M_PI/180.; 34 m_RadLa=latitude*M_PI/180.; 35 Ec=RJ+(RC-RJ)*(90.-m_Latitude)/90.; 36 Ed=Ec*cos(m_RadLa); 37 } 38 39 return self; 40 } 41 @end
1 #import "ViewController.h" 2 #import "MyLatLng.h" 3 4 @interface ViewController () 5 6 @end 7 8 @implementation ViewController 9 10 - (void)viewDidLoad { 11 [super viewDidLoad]; 12 // Do any additional setup after loading the view, typically from a nib. 13 14 MyLatLng *A = [[MyLatLng alloc]init:113.249648 latitude:23.401553]; 15 MyLatLng *B = [[MyLatLng alloc]init:113.246033 latitude:23.403362]; 16 17 [self getAngle:A B:B]; 18 } 19 20 - (double)getAngle:(MyLatLng *)A B:(MyLatLng *)B { 21 double dx=(B.m_RadLo-A.m_RadLo)*A.Ed; 22 double dy=(B.m_RadLa-A.m_RadLa)*A.Ec; 23 double angle=0.0; 24 angle = atan(fabs(dx/dy))*180./M_PI; 25 double dLo=B.m_Longitude-A.m_Longitude; 26 double dLa=B.m_Latitude-A.m_Latitude; 27 if(dLo>0&&dLa<=0){ 28 angle=(90.-angle)+90; 29 } 30 else if(dLo<=0&&dLa<0){ 31 angle=angle+180.; 32 }else if(dLo<0&&dLa>=0){ 33 angle= (90.-angle)+270; 34 } 35 return angle; 36 } 37 38 - (void)didReceiveMemoryWarning { 39 [super didReceiveMemoryWarning]; 40 // Dispose of any resources that can be recreated. 41 }