clear;clc;
addpath /home/yanjiang/matbin;
cd /home/yanjiang/gstatdir;
load alldata.mat;
cd /chan-data1/yjiang/SIM;
[lont latt vt sigt]=textread('180_truth','%n%n%n%n');
x=alldata(:,1);y=alldata(:,2);v=alldata(:,3);sigv=alldata(:,4);
[dlat,dlon]=lg2dg(y.*1000,x.*1000,deg2rad(75),0);
lat=rad2deg(dlat+deg2rad(75));
lon=rad2deg(dlon+deg2rad(320));
alldata2=[lon lat v sigv];pdmt=[];
for m=round(min(lon)):1:round(max(lon))
    for n=round(min(lat)):1:round(max(lat))
        alldata3=alldata2;lon=alldata3(:,1);
        alldata3(any(lon<m-0.5 | lon>m+0.5,2),:)=[];
        lat=alldata3(:,2);
        alldata3(any(lat<n-0.5 | lat>n+0.5,2),:)=[];
        v3=alldata3(:,3);sigv3=alldata3(:,3);
        %av=mean(v3);
        av=sum(v3.*sigv3.^(-2))./sum(sigv3.^(-2));
        sigav=1./sum(sigv3.^(-2));
        truth=[lont latt vt sigt];
        truth(any(lont~=m | latt~=n,2),:)=[];
        tv=truth(:,3);
        dv=abs((av-tv)./tv).*100;
        dv2=abs(av-tv);
        pdmt1=[m n av sigav dv dv2 tv];
        pdmt=[pdmt;pdmt1];
    end
end
pdmt(any(pdmt(:,5)>200,2),:)=[];pdmt(any(isnan(pdmt(:,5)),2),:)=[];
lon=pdmt(:,1);lat=pdmt(:,2);v=pdmt(:,3);sigv=pdmt(:,4);dv2=pdmt(:,6);
figure (1); scatter(lon,lat,25,dv2,'.');colormap (jet);hold on;colorbar;
figure (2); scatter(lon,lat,25,pdmt(:,7),'.');
colormap (jet);hold on;colorbar;
truth=[lont latt vt sigt];
for i=1:length(pdmt)
    truth(any(truth(:,1)-pdmt(i,1)==0 & truth(:,2)-pdmt(i,2)==0,2),:)=[];
end
pdmttruth=[pdmt(:,1) pdmt(:,2) pdmt(:,3) pdmt(:,4)];
truth=[truth;pdmttruth];
    
    
    
    
    
    
    
    
    
    
    
    
    