lat_MODIS=-89.5:1:89.5;
lon_MODIS=-179.5:1:179.5;
load /Users/yang146/Aerosol_Variaiton/IRUN/AODBD/ncfile/corr_AODBD_MODIS_WINDGobi.mat 
load ~/NC/grid.mat
load /Users/yangyang/Aerosol_Variaiton/IRUN/AODBD/ncfile/AOD_CEDS.mat

AOD_ANTH=mean(AODBC_DJF(:,:,2:6)+AODPOM_DJF(:,:,2:6)+AODSOA_DJF(:,:,2:6)+AODSO4_DJF(:,:,2:6),3);
AOD_DUST=mean(AODDUST_DJF(:,:,2:6),3);

[~,lonmin]=min(abs(lon-100));
[~,lonmax]=min(abs(lon-110));
[~,latmin]=min(abs(lat-35));
[~,latmax]=min(abs(lat-45));

n=0;
for i=lonmin:lonmax
    for j=latmin:latmax
        n=n+1;
%        F_ANTH(n)=AOD_ANTH(i,j)/(AOD_ANTH(i,j)+AOD_DUST(i,j));
%        F_DUST(n)=AOD_DUST(i,j)/(AOD_ANTH(i,j)+AOD_DUST(i,j));
        
        F_ANTH(n)=AOD_ANTH(i,j);
        F_DUST(n)=AOD_DUST(i,j);
        
        [~,lon_MODIS2]=min(abs(lon_MODIS-lon(i)));
        [~,lat_MODIS2]=min(abs(lat_MODIS-lat(j)));
        
        R_MODIS(n)=R(lon_MODIS2,lat_MODIS2);
    end
end
        

X = [ones(size(F_ANTH')) F_ANTH' F_DUST'];
y = R_MODIS';

b = regress(y,X);