function OK = convertCoordinates_test_netherlands(varargin)
%CONVERTCOORDINATES_TEST_NETHERLANDS  test for Dutch points
%
%    OK = convertcoordinates_test_netherlands(<EPSG>)
%
%  x   is accurate to m
%  y   is accurate to m
%  lon is accurate to deg
%  lat is accurate to deg
%
%See also: convert_coordinates

   if nargin==1
   EPSG = varargin{1};
   else
   EPSG = load('EPSG');
   end

   OK(1)       = test3(EPSG);
   OK(2)       = test4(EPSG);
   OK          = all(OK);

function OK = test3(EPSG)
%  Van een Kernnetpunt
%  https://rdinfo.kadaster.nl/?inhoud=/rd/info.html%23publicatie&navig=/rd/nav_serverside.html%3Fscript%3D1
%  https://rdinfo.kadaster.nl/pics/publijst1.gif

D.Puntnummer        = '019111';
D.Actualiteitsdatum = datenum(1999,6,1);
D.Nr                = 17;
D.X                 = 155897.26;
D.Y                 = 603783.39;
D.H                 = 3.7;
D.NB                = 53+25/60+13.2124/3600;
D.OL                = 05+24/60+02.5391/3600;
D.h                 = 44.83;

[lon,lat] = convertCoordinates(D.X ,D.Y ,EPSG,'CS1.code',28992,'CS2.code', 4326);
[X  ,Y  ] = convertCoordinates(D.OL,D.NB,EPSG,'CS1.code', 4326,'CS2.code',28992);
[X2 ,Y2 ] = convertCoordinates(lon ,lat ,EPSG,'CS1.code', 4326,'CS2.code',28992); % and back

% WGS84 and ETRS89 are not identical. WGS84 is < 1 m accurate
% The difference in 2004 is say 35 centimeter, see http://www.rdnap.nl/stelsels/stelsels.html
% So for testing less < 0.5 m error is OK.

% num2str(D.OL - lon) check projection onesided
% num2str(D.NB - lat)
%
% num2str(D.X - X)    check projection onesided
% num2str(D.Y - Y)
%
% num2str(D.X - X2)   check projection twosided: internal consistensy
% num2str(D.Y - Y2)

OK = abs(D.OL - lon) < 1e-5 & abs(D.NB - lat) < 1e-5 & ...
    abs(X -D.X)     < 0.5  & abs(Y -D.Y)     < 0.5  & ...
    abs(X2-D.X)     < 0.5  & abs(Y2-D.Y)     < 0.5;

function OK = test4(EPSG)
%% comparison with RDNAPTRANS reference coordinates

% Errors should be smaller than maximum the distortion in the RD system,
% which is less than 25 cm.  
% To convert coordinates including the distortion in the RD grid a
% dedicated program with the distortion grid must be used, such as
% RDNAPTRANS

%from: http://www.kadaster.nl/index_frames.html?inhoud=/rijksdriehoeksmeting/homepage.html&navig=/rijksdriehoeksmeting/nav_serverside.html%3Fscript%3D1
% Appendix		RDNAPTRANSTM2008 test sheet 
% 
% The RDNAPTRANSTM2008 procedure should be tested in both directions. 
% The differences with the coordinates listed below should not exceed:
% RD x and y coordinates:                      0.001      meters
% NAP heights and ETRS89 ellipsoidal heights:  0.001      meters
% ETRS89 latitude and longitude:               0.00000001 degrees
% 
% From ETRS89 to RD/NAP
%                         |-----------ETRS89--------------------|--------------RD/NAP----------------|                 		
% No.   Name              | latitude (°)  longitude (°)  h (m)  |      x (m)        y (m)    NAP (m) |
% 01    Texel             | 53.160753042  4.824761912   42.8614 | 117380.1200  575040.3400    1.0000 |
% 02    Noord-Groningen   | 53.419482050  6.776726674   42.3586 | 247380.5600  604580.7800    2.0000 |
% 03    Amersfoort        | 52.155172897  5.387203657   43.2551 | 155000.0000  463000.0000    0.0000 |
% 04    Amersfoort 100m   | 52.155172910  5.387203658  143.2551 | 155000.0000  463000.0000  100.0000 |
% 05    Zeeuws-Vlaanderen | 51.368607152  3.397588595   47.4024 |  16460.9100  377380.2300    3.0000 |
% 06    Zuid-Limburg      | 50.792584916  5.773795548  245.9478 | 182260.4500  311480.6700  200.0000 |
% 07    Maasvlakte        | 51.947393898  4.072887101   47.5968 |  64640.8900  440700.0101    4.0000 |
% 08*   outside           | 48.843030210  8.723260235   52.0289 | 400000.2300  100000.4500    5.0000 |
% 09*   no_rd&geoid       | 50.687420392  4.608971813   51.6108 | 100000.6700  300000.8900    6.0000 |
% 10*   no_geoid          | 51.136825197  4.601375361   50.9672 | 100000.6700  350000.8900    6.0000 |
% 11*   no_rd             | 52.482440839  4.268403889   49.9436 |  79000.0100  500000.2300    7.0000 |
% 12*   edge_rd           | 51.003976532  3.891247830   52.7427 |  50000.4500  335999.6700    8.0000 |
% 

ETRS89toRDNAP = [
    53.160753042  4.824761912   42.8614  117380.1200  575040.3400    1.0000
    53.419482050  6.776726674   42.3586  247380.5600  604580.7800    2.0000
    52.155172897  5.387203657   43.2551  155000.0000  463000.0000    0.0000
    52.155172910  5.387203658  143.2551  155000.0000  463000.0000  100.0000
    51.368607152  3.397588595   47.4024   16460.9100  377380.2300    3.0000
    50.792584916  5.773795548  245.9478  182260.4500  311480.6700  200.0000
    51.947393898  4.072887101   47.5968   64640.8900  440700.0101    4.0000
    48.843030210  8.723260235   52.0289  400000.2300  100000.4500    5.0000
    50.687420392  4.608971813   51.6108  100000.6700  300000.8900    6.0000
    51.136825197  4.601375361   50.9672  100000.6700  350000.8900    6.0000
    52.482440839  4.268403889   49.9436   79000.0100  500000.2300    7.0000
    51.003976532  3.891247830   52.7427   50000.4500  335999.6700    8.0000
    ];


%% from ETRS 89 to RD

    [RDx,RDy,OPT] = convertCoordinates(ETRS89toRDNAP(1,2),ETRS89toRDNAP(1,1),EPSG,...
        'CS1.name', 'ETRS89','CS1.type','geographic 2D','CS2.code',28992);
    % from the OPT, several alternative methods for the datum transformation
    % can be retrieved. Try them all
    max_error = [];
    for ii = 1:length(OPT.datum_trans.alt_code );
        
        [RDx,RDy] = convertCoordinates(ETRS89toRDNAP(:,2),ETRS89toRDNAP(:,1),EPSG,...
            'CS1.name', 'ETRS89','CS1.type','geographic 2D','CS2.code',28992,...
            'datum_trans.code',OPT.datum_trans.alt_code(ii));
        max_error(end+1)  =  max(((RDx - ETRS89toRDNAP(:,4)).^2+(RDy - ETRS89toRDNAP(:,5)).^2).^.5);
        % mean_error = mean(((RDx - ETRS89toRDNAP(:,4)).^2+(RDy - ETRS89toRDNAP(:,5)).^2).^.5);

        % fprintf(1,'ETRS89 to RD using method % 25s, max error =% 8.4f m, mean error =% 8.4f m. Deprecated = %s\n',...
        %    OPT.datum_trans.alt_name{ii},max_error,mean_error,OPT.datum_trans.alt_deprecated{ii});
    end
    
    % fprintf('ETRS89 to RD default: % 29s\n\n',OPT.datum_trans.name{1})
    
    % OK is if all conversions give a max error smaller than 25cm
    OK(1) = all(max_error < 0.25);
%% from WGS 84 to RD

    [RDx,RDy,OPT] = convertCoordinates(ETRS89toRDNAP(1,2),ETRS89toRDNAP(1,1),EPSG,...
        'CS1.name', 'WGS 84','CS1.type','geographic 2D','CS2.code',28992);
    % from the OPT, several alternative methods for the datum transformation
    % can be retrieved. Try them all
    max_error = [];
    for ii = 1:length(OPT.datum_trans.alt_code );
        [RDx,RDy] = convertCoordinates(ETRS89toRDNAP(:,2),ETRS89toRDNAP(:,1),EPSG,...
            'CS1.name', 'WGS 84','CS1.type','geographic 2D','CS2.code',28992,...
            'datum_trans.code',OPT.datum_trans.alt_code(ii));
        max_error(end+1)  =  max(((RDx - ETRS89toRDNAP(:,4)).^2+(RDy - ETRS89toRDNAP(:,5)).^2).^.5);
%         mean_error = mean(((RDx - ETRS89toRDNAP(:,4)).^2+(RDy - ETRS89toRDNAP(:,5)).^2).^.5);
%         
%         fprintf(1,'WGS 84 to RD using method % 25s, max error =% 8.4f m, mean error =% 8.4f m. Deprecated = %s\n',...
%             OPT.datum_trans.alt_name{ii},max_error,mean_error,OPT.datum_trans.alt_deprecated{ii});
    end
%     fprintf('WGS 84 to RD default: % 29s\n\n',OPT.datum_trans.name{1})

% OK is if all but the first conversion give a max error smaller than 25cm
OK(2) = all(max_error(2:end) < 0.25);

%% From RD/NAP to ETRS89
% 
% No.    Name    RD/NAP            ETRS89
%                         |-----------ETRS89-------------------|--------------RD/NAP------------------|
%                         |      x (m)        y (m)    NAP (m) | latitude (°)  longitude (°)   h (m)  |
% 01    Texel             | 117380.1200  575040.3400    1.0000 | 53.160753042   4.824761912   42.8614 |
% 02    Noord-Groningen   | 247380.5600  604580.7800    2.0000 | 53.419482050   6.776726674   42.3586 |
% 03    Amersfoort        | 155000.0000  463000.0000    0.0000 | 52.155172897   5.387203657   43.2551 |
% 04    Amersfoort_100m   | 155000.0000  463000.0000  100.0000 | 52.155172910   5.387203658  143.2551 |
% 05    Zeeuws-Vlaanderen |  16460.9100  377380.2300    3.0000 | 51.368607152   3.397588595   47.4024 |
% 06    Zuid-Limburg      | 182260.4500  311480.6700  200.0000 | 50.792584916   5.773795548  245.9478 |
% 07    Maasvlakte        |  64640.8900  440700.0100    4.0000 | 51.947393898   4.072887101   47.5968 |
% 08*   outside           | 400000.2300  100000.4500    5.0000 | 48.843030210   8.723260235   52.0289 |
% 09*   no_rd&geoid       | 100000.6700  300000.8900    6.0000 | 50.687420392   4.608971813   51.6108 |
% 10*   no_geoid          | 100000.6700  350000.8900    6.0000 | 51.136825197   4.601375361   50.9672 |
% 11*   no_rd             |  79000.0100  500000.2300    7.0000 | 52.482440839   4.268403889   49.9436 |
% 12*   edge_rd           |  50000.4500  335999.6700    8.0000 | 51.003976532   3.891247830   52.7427 |
% 							      
% 
% *) Points 08 - 12 are outside the region where interpolation between either the NLGEO2004 geoid or the RD correction grid points is possible. If coordinates are computed for these points, the output should be accompanied by a warning.


RDNAPtoETRS89 = [
	117380.1200  575040.3400    1.0000  53.160753042   4.824761912   42.8614
	247380.5600  604580.7800    2.0000  53.419482050   6.776726674   42.3586
	155000.0000  463000.0000    0.0000  52.155172897   5.387203657   43.2551
	155000.0000  463000.0000  100.0000  52.155172910   5.387203658  143.2551
	 16460.9100  377380.2300    3.0000  51.368607152   3.397588595   47.4024
	182260.4500  311480.6700  200.0000  50.792584916   5.773795548  245.9478
	 64640.8900  440700.0100    4.0000  51.947393898   4.072887101   47.5968
	400000.2300  100000.4500    5.0000  48.843030210   8.723260235   52.0289
	100000.6700  300000.8900    6.0000  50.687420392   4.608971813   51.6108
	100000.6700  350000.8900    6.0000  51.136825197   4.601375361   50.9672
	 79000.0100  500000.2300    7.0000  52.482440839   4.268403889   49.9436
	 50000.4500  335999.6700    8.0000  51.003976532   3.891247830   52.7427
	 ];

 %% from RD to ETRS89

     [lon,lat,OPT] = convertCoordinates(RDNAPtoETRS89(1,1),RDNAPtoETRS89(1,2),EPSG,...
         'CS2.name', 'ETRS89','CS2.type','geographic 2D','CS1.code',28992);
     % from the OPT, several alternative methods for the datum transformation
     % can be retrieved. Try them all
     max_error = [];
     for ii = 1:length(OPT.datum_trans.alt_code );
         [lon,lat] = convertCoordinates(RDNAPtoETRS89(:,1),RDNAPtoETRS89(:,2),EPSG,...
             'CS2.name', 'ETRS89','CS2.type','geographic 2D','CS1.code',28992,...
             'datum_trans.code',OPT.datum_trans.alt_code(ii));
         max_error(end+1)  =  max(((lat - RDNAPtoETRS89(:,4)).^2+(lon - RDNAPtoETRS89(:,5)).^2).^.5);
%          mean_error = mean(((lat - RDNAPtoETRS89(:,4)).^2+(lon - RDNAPtoETRS89(:,5)).^2).^.5);
         
%          fprintf(1,'RD to ETRS89 using method % 25s, max error =% 10.8f deg, mean error =% 10.8f deg. Deprecated = %s\n',...
%              OPT.datum_trans.alt_name{ii},max_error,mean_error,OPT.datum_trans.alt_deprecated{ii});
     end
%      fprintf('RD to ETRS89 default: % 29s\n\n',OPT.datum_trans.name{1})

OK(3) = all(max_error < 1e-5);

 %% from RD to WGS 84

     [lon,lat,OPT] = convertCoordinates(RDNAPtoETRS89(1,1),RDNAPtoETRS89(1,2),EPSG,...
         'CS2.name', 'WGS 84','CS2.type','geographic 2D','CS1.code',28992);
     % from the OPT, several alternative methods for the datum transformation
     % can be retrieved. Try them all
     max_error = [];
     for ii = 1:length(OPT.datum_trans.alt_code );
         [lon,lat] = convertCoordinates(RDNAPtoETRS89(:,1),RDNAPtoETRS89(:,2),EPSG,...
             'CS2.name', 'WGS 84','CS2.type','geographic 2D','CS1.code',28992,...
             'datum_trans.code',OPT.datum_trans.alt_code(ii));
         max_error(end+1)  =  max(((lat - RDNAPtoETRS89(:,4)).^2+(lon - RDNAPtoETRS89(:,5)).^2).^.5);
%          mean_error = mean(((lat - RDNAPtoETRS89(:,4)).^2+(lon - RDNAPtoETRS89(:,5)).^2).^.5);
%          
%          fprintf(1,'RD to WGS 84 using method % 25s, max error =% 10.8f deg, mean error =% 10.8f deg. Deprecated = %s\n',...
%              OPT.datum_trans.alt_name{ii},max_error,mean_error,OPT.datum_trans.alt_deprecated{ii});
     end
%      fprintf('RD to WGS 84 default: % 29s\n\n',OPT.datum_trans.name{1})
% OK is if all but the first conversion give a max error smaller 1e-5
OK(4) = all(max_error(2:end) < 1e-5);

%% 
OK = all(OK);
