I tried to calculate the velocity with linear regression from GPS data with the following script but there is an error

%Titik Ikat
FID = fopen('FIXSTA.SIG');
datfix = textscan(FID,'%s');

Error using textscan
Invalid file identifier. Use fopen to generate a valid file identifier.

fclose(FID);
tolsig=(size(datfix{1,1},1)-15)/4;
for ts=1:tolsig
    FCS{ts,1}=datfix{1,1}{(ts-1)*4+16,1};
end

%Poligon defined area
polg=[batss(1,1) batss(3,1); batss(1,1) batss(4,1);batss(1,1) batss(2,1); batss(2,1) batss(4,1)];


cd(dirasr);

%IMPORT Big Data
FID = fopen(nafilout);
dataz = textscan(FID,'%s');
fclose(FID);
redat=reshape(dataz{1,1},9,[])';
natik=redat(:,1);
jei=1; awa=0;
wakt=0;
fff = waitbar(0,['Reading Raw Files']);
for qq1=1:size(redat,1)
    tic
    for qq2=2:size(redat,2)
        comem(qq1,qq2-1)=str2num(char(redat{qq1,qq2}));
    end
    if qq1<size(redat,1)
        if strcmpi(natik(qq1,1),natik(qq1+1,1))
            continue;
        else
            %SPLIT Data
            ntk(jei,1)=natik(qq1,1);
            dat{jei}=comem(awa+1:qq1,:);
            awa=qq1;
            jei=jei+1;
        end
    end     
    
    wakt=wakt+toc;
    perc=qq1/(size(redat,1)-1+1)*100; 
    if wakt/60<1;  waktu=[num2str(wakt,'%.3f'),' sec'];
    else
        fwk=floor(wakt/60);
        if wakt/3600<1;    waktu=[num2str(fwk,'%.3f'),' min ',num2str(wakt-fwk*60,'%.3f'),' sec'];
        else waktu=[num2str(floor(fwk/60),'%.3f'),' hour ', num2str(fwk-floor(fwk/60)*60,'%.3f'),' min ',num2str(wakt-fwk*60,'%.3f'),' sec'];
        end    
    end    
    waitbar(perc/100,fff,['Reading Raw Files ', num2str(perc,'%.3f'), ' % (',waktu,')']);
end
close(fff)

ntk(jei,1)=natik(end,1);
dat{jei}=comem(awa+1:end,:);
daseb=size(ntk,1);

cd(diras);

%File Offset
if isfile('offset.txt')==1
    FID = fopen(nafiloff);
    dataoff = textscan(FID,'%s');
    fclose(FID);
    matoff=reshape(dataoff{1,1},3,[])';
else
    matoff=0;
end

ros=1;

spc=[' '];
start=1;
endd=daseb;

wakt=0;
fff = waitbar(0,['Processing sites']);
for j=start:endd
    tic
    SPCE(j,:)=spc;
    
    
    %Offset per file
    clear funcch funcof
    matoffeach=matoff(find(strcmpi(matoff(:,1),ntk{j})),:);
    if isempty(matoffeach)==0
        mfo1=find(strcmpi(matoffeach(:,3),'A'));
        mfo2=find(strcmpi(matoffeach(:,3),'B'));
        mfo3=find(strcmpi(matoffeach(:,3),'C'));
        funcof=matoffeach([mfo1;mfo3],:);
        funcch=matoffeach([mfo2;mfo3],:);
        if isempty(funcof)==0
            indiko=1; 
        else
            indiko=0;
        end
        if isempty(funcch)==0
            indikc=1; 
        else
            indikc=0;
        end
        
    else
        indiko=0;
        indikc=0;
    end
    
    %Titik Ikat
    SFCS=sum(ismember(ntk{j},FCS));
    if SFCS==0
        igsgps=0;
    else
        igsgps=1;
    end
    
    
    %Pendefinisian koordinat kartesian geosentrik
    xg=dat{j}(:,3);
    yg=dat{j}(:,4);
    zg=dat{j}(:,5);
    sd=size(xg,1);          %Jumlah data per titik
    
    for kckc=1:sd
        %Pendefinisian koordinat geodetik
        [lg,bg,hg]=cg2geod(xg(kckc,1),yg(kckc,1),zg(kckc,1));
    end
    
    xrg=mean(xg);
    yrg=mean(yg);
    zrg=mean(zg);
    
    dxg=xg-xrg*ones(sd,1);
    dyg=yg-yrg*ones(sd,1);
    dzg=zg-zrg*ones(sd,1);
    
    lrg=mean(lg);
    brg=mean(bg);
    lin(j,1)=lrg;
    buj(j,1)=brg; 
    
    %Inside area
    statinsp=inpolygon(brg,lrg,polg(:,1),polg(:,2));
    
    %Pendefinisian koordinat kartesian toposentrik
    ROTASI=[-sind(lrg)*cosd(brg) -sind(lrg)*sind(brg)   cosd(lrg)
        -sind(brg)         cosd(brg)              0
        cosd(lrg)*cosd(brg) cosd(lrg)*sind(brg)   sind(lrg)];
    for k=1:sd
        neug=ROTASI*[dxg(k,1);dyg(k,1);dzg(k,1)];
        ng(k,1)=neug(1,1);
        eg(k,1)=neug(2,1);
        ug(k,1)=neug(3,1);
    end  
    neu=[ng eg ug]*1000;    %m > mm
    
    %Pendefinisian ketelitian toposentrik
    ketn=dat{j}(:,6)+0.0000001;
    kete=dat{j}(:,7)+0.0000001;
    ketu=dat{j}(:,8)+0.0000001;
    
    %Hapus outlier dari ketelitian
    fk1=find(ketn>0.01);
    fk2=find(kete>0.01);
    fk3=find(ketu>0.01);
    fk=unique([fk1;fk2;fk3]);
    ketn(fk,:)=[];
    kete(fk,:)=[];
    ketu(fk,:)=[];
    neu(fk,:)=[];
    ket=[ketn kete ketu]*1000;  %m > mm
    
    sd=size(neu,1);
    
    %Pendefinisian sumbu waktu
    tah=dat{j}(:,1)+2000;
    doy=dat{j}(:,2);
    sumx=tah+doy/365;
    sumx(fk,:)=[];
    
    if bamixall==1
        bamin=min(sumx)-(max(sumx)-min(sumx))/20;
        bamax=max(sumx)+(max(sumx)-min(sumx))/20;
    else
        bamin=bam1;
        bamax=bam2;
        fint1=find(sumx>bamin);
        fint2=find(sumx<bamax);
        fint3=intersect(fint1,fint2);
        sumx=sumx(fint3,1);
        if isempty(sumx)==1
            continue;
        end
        ketn=ketn(fint3,:);
        kete=kete(fint3,:);
        ketu=ketu(fint3,:);
        neu=neu(fint3,:);
        sd=size(neu,1);
    end
    
    %Pendefinisian matrik yang generic untuk perataan parameter
    matB=[ones(sd,1) sumx];
    
    %STEP FUNCTION
    sks=1;
    if indiko==1
        sfcf=size(funcof,1);
        for inST=1:sfcf
            doon(sks,1)=str2num(funcof{inST,2});
            sks=sks+1;
            matB(:,inST+2)=ones(sd,1);
            matB(matB(:,2)<str2num(funcof{inST,2}),inST+2)=0;
        end
    end
    if indikc==1
        if indiko==0;
            sfcf=0;
        end
        sks2=sfcf+size(funcch,1);
        for inST2=1:size(funcch,1)
            doon(sks,1)=str2num(funcch{inST2,2});
            sks=sks+1;
            matB(:,inST2+2+sfcf)=sumx-str2num(funcch{inST2,2});
            matB(matB(:,2)<str2num(funcch{inST2,2}),inST2+2+sfcf)=0;
        end
    else
        if indiko==1
            sks2=sfcf;
        else    sks2=0;     end
    end
    
    SNN=1./ketn.^2;
    SNE=1./kete.^2;
    SNU=1./ketu.^2;
    Z=[SNN SNE SNU];
    
    wakt=wakt+toc;
    perc=(j-0.7)/(endd-start+1)*100;
    if wakt/60<1;  waktu=[num2str(wakt,'%.3f'),' sec'];
    else
        fwk=floor(wakt/60);
        if wakt/3600<1;    waktu=[num2str(fwk,'%.3f'),' min ',num2str(wakt-fwk*60,'%.3f'),' sec'];
        else waktu=[num2str(floor(fwk/60),'%.3f'),' hour ', num2str(fwk-floor(fwk/60)*60,'%.3f'),' min ',num2str(wakt-fwk*60,'%.3f'),' sec'];
        end    
    end    
    waitbar(perc/100,fff,['Processing sites ', num2str(perc,'%.3f'), ' % (',waktu,')']);
    tic
    
    %Penetuan jumlah komponen (NS/EW/UD)
    nuof=3;
    
    %Perhitungan per komponen
    f=figure('visible','off');
    thresd=2*std(neu);
    for i=1:nuof
        
        %Pendefinisian matrik yang spesifik untuk perataan parameter
        F=-[neu(:,i)];
        P=diag(Z(:,i));
        
        %Perhitungan hitung perataan parameter
        qxx=inv(matB'*P*matB);
        akhir=-(qxx)*matB'*P*F;
        VV=matB*akhir+F;
        varbak=(VV'*P*VV)/(sd-2);
        varpar= varbak*qxx;
        
        %Remove outlier
        thres=1*std(detrend(neu(:,i)));
        neu(abs(VV(:,1))>(thres),:)=[];
        F(abs(VV(:,1))>(thres),:)=[];
        sumx(abs(VV(:,1))>(thres),:)=[];
        ket(abs(VV(:,1))>(thres),:)=[];
        matB(abs(VV(:,1))>(thres),:)=[];
        Z(abs(VV(:,1))>(thres),:)=[];
        P(abs(VV(:,1))>(thres),:)=[];
        P(:,abs(VV(:,1))>(thres))=[];
        sd2=size(neu,1);
        
    end
    
    msz=(floor(2/log(nthroot(size(matB,1),6)))+1)^2;
    if msz>10
        tipe=['.'];
        msz=20;
    else
        tipe=['o'];
        msz=2;
    end
    
    
    wakt=wakt+toc;
    perc=(j-0.6)/(endd-start+1)*100;
    if wakt/60<1;  waktu=[num2str(wakt,'%.3f'),' sec'];
    else
        fwk=floor(wakt/60);
        if wakt/3600<1;    waktu=[num2str(fwk,'%.3f'),' min ',num2str(wakt-fwk*60,'%.3f'),' sec'];
        else waktu=[num2str(floor(fwk/60),'%.3f'),' hour ', num2str(fwk-floor(fwk/60)*60,'%.3f'),' min ',num2str(wakt-fwk*60,'%.3f'),' sec'];
        end    
    end    
    waitbar(perc/100,fff,['Processing sites ', num2str(perc,'%.3f'), ' % (',waktu,')']);
    tic
    
    for i=1:nuof
        %Pendefinisian matrik yang spesifik untuk perataan parameter
        F=-[neu(:,i)];
        P=diag(Z(:,i));
        
        %Perhitungan hitung perataan parameter
        akhir=-(inv(matB'*P*matB))*matB'*P*F;
        VV=matB*akhir+F;
        varbak=(VV'*P*VV)/(sd2-2);
        qxx=inv(matB'*P*matB);
        varpar= varbak*qxx;
        
        %Hasil akhir berupa kecepatan dan ketelitiannya
        AK1(nuof+1-i)=akhir(2,1);
        AK2(nuof+1-i)=sqrt(varpar(2,2));
        
        if exist("sfcf")==1
            AKO(nuof+1-i,:)=[akhir(3:sfcf+2,1)]';
            perako=sqrt(varpar(3:3:sfcf+2,3:3:sfcf+2));
            AKO2(nuof+1-i,:)=diag(perako);
            MAKO(nuof+1-i)=max(AKO(nuof+1-i,:));
            MAKO2(nuof+1-i)=max(AKO2(nuof+1-i,:));
            if sks2>0
                for itam=1:sks2-sfcf
                    AKP{itam}(nuof+1-i)=akhir(itam+sfcf+2,1);
                end
            end
        end
        
        %Pembuatan garis regresi linear
        ssux=[bamin;sumx;bamax];
        if exist("sks2")==1
            if sks2>0
                for k=1:size(doon,1)
                    ssux=[ssux;doon(k,1);doon(k,1)+0.001];
                end
                ssux=sortrows(ssux,1);
            end
        end
        sssux=[bamin;bamax];
        
        sumxx{i}=sumx;
        yperb=akhir(1,1)+akhir(2,1)*sumx;
        yperb2=yperb;
        y=akhir(1,1)+akhir(2,1)*ssux;
        if indiko==1
            for inST=1:sfcf
                y(ssux(:,1)>str2num(funcof{inST,2}),:)=y(ssux(:,1)>str2num(funcof{inST,2}),:)+akhir(inST+2,1);
                yperb(sumx(:,1)>str2num(funcof{inST,2}),:)=yperb(sumx(:,1)>str2num(funcof{inST,2}),:)+akhir(inST+2,1);
            end
        end
        if indikc==1
            for inST2=1:size(funcch,1)
                y(ssux(:,1)>str2num(funcch{inST2,2}),:)=y(ssux(:,1)>str2num(funcch{inST2,2}),:)+akhir(inST2+2+sfcf,1)*(ssux(ssux(:,1)>str2num(funcch{inST2,2}),1)-str2num(funcch{inST2,2}));
                yperb(sumx(:,1)>str2num(funcch{inST2,2}),:)=yperb(sumx(:,1)>str2num(funcch{inST2,2}),:)+akhir(inST2+2+sfcf,1)*(sumx(sumx(:,1)>str2num(funcch{inST2,2}),1)-str2num(funcch{inST2,2}));
            end
        end
        
        yy=[y(1,1);y(end,1)];
        resd{i}=neu(:,i)-yperb;
        %             neu(:,i)=yperb+resd{i}/2;
        %             resd{i}=neu(:,i)-yperb;
        
        %Subplot tiga baris untuk metode 1
        subplot(3,1,i)
        POS=get(subplot(3,1,i), 'Position' );
        POS(2) =POS(2)+0.05*(i-2);                       % Position move 
        set(subplot(3,1,i), 'Position', POS);
        hold on
        grid on
        
        %Plot ketelitian
        for rr=1:sd2
            plot([sumx(rr,1) sumx(rr,1)],[neu(rr,i)-ket(rr,i) neu(rr,i)+ket(rr,i)],'green','LineWidth',0.2)
        end
        
        %Plot garis regresi linear dan titik koordinat
        plot(ssux,y,'--black','LineWidth',1);
        plot(sumx,neu(:,i),tipe,'color','red','MarkerSize',msz);
        
        %Plot label. dan batas sumbu
        xlim([bamin bamax]);
        if fixylim==0
            ylim([minp1 minp2]);
            iminp=(minp2-minp1)/6;
            set(gca,'YTick',[minp1:iminp:minp2]);
        else
            ylim([min(neu(:,i))-1 max(neu(:,i))+1])
        end
        if i==1
            title(['Site: ',char(ntk(j,:))],'FontSize',14);
            ylabel('North (mm)','FontSize',14);
            set(gca,'XTickLabel',[])
        elseif i==2
            ylabel('East (mm)','FontSize',14);
            set(gca,'XTickLabel',[])
        elseif i==3
            xlabel('Year','FontSize',14);
            ylabel('Up (mm)','FontSize',14);
        end
    end
    %SELESAI tiga komponen per titik
    
    %save neu
    neudatb{j}=[sumx neu sqrt(1./Z)];
    neu4save=neudatb{j};
    nafil=[num2str(j,'%03.f') '-' ntk{j} '.txt'];
    cd(dirts)
    if isdir('Coor-Topo')==0
        mkdir Coor-Topo
    end
    cd Coor-Topo
    save(nafil,'neu4save','-ascii');
    cd ..
    %---------------------------------------------------------------------------------------------------------------------
    
    %Penyimpanan figure
    cd(dirts)
    if isdir('TS')==0
        mkdir TS
    end
    cd TS
    saveas(gcf,['Time Series of ',char(ntk(j,:)),'.png'])
    cd ..
    
    f=figure('visible','off');
    for gg=1:nuof
        subplot(3,1,gg)
        
        POS=get(subplot(3,1,gg), 'Position' );
        POS(2) =POS(2)+0.05*(gg-2);                       % Position move 
        set(subplot(3,1,gg), 'Position', POS);
        
        hold on
        grid on
        
        %Plot garis regresi linear dan titik koordinat
        plot([bamin bamax],[0 0],'--black','LineWidth',1);
        plot(sumxx{gg},resd{gg},tipe,'color','red','MarkerSize',msz);
        JRES{j,gg}=resd{gg};
        JSUM{j,gg}=sumxx{gg};
        
        %Plot label. dan batas sumbu
        xlim([bamin bamax]);
        if fixylim==0
            ylim([minr1 minr2]);
            iminr=(minr2-minr1)/4;
            set(gca,'YTick',[minr1:iminr:minr2]);
        else
            ylim([min(resd{gg})-1 max(resd{gg})+1])
        end
        if gg==1
            title(['Site: ',char(ntk(j,:))],'FontSize',14);
            ylabel('North (mm)','FontSize',14);
            set(gca,'XTickLabel',[])
        elseif gg==2
            ylabel('East (mm)','FontSize',14);
            set(gca,'XTickLabel',[])
        elseif gg==3
            xlabel('Year','FontSize',14);
            ylabel('Up (mm)','FontSize',14);
        end
        
    end
    if isdir('RTS')==0
        mkdir RTS
    end
    cd RTS
    saveas(gcf,['Residual Time Series of ',char(ntk(j,:)),'.png'])
    cd ..
    
    f=figure('visible','off');
    
    cd(diras)
    
    
    %Penyimpanan data untuk plot GMT
    OOOO(j,1:2)=[brg lrg];
    if jenp==1
        OOOO(j,3:7)=[AK1(:,1+nuof-2:2+nuof-2) AK2(:,1+nuof-2:2+nuof-2) 0];
        VOO(j,:)=[OOOO(j,1:2) 0 AK1(:,1) 0 AK2(:,1) 0];
    elseif jenp==2
        OOOO(j,3:7)=[MAKO1(:,1+nuof-2:2+nuof-2) MAKO2(:,1+nuof-2:2+nuof-2) 0];
        VOO(j,:)=[OOOO(j,1:2) 0 MAKO(:,1) 0 MAKO2(:,1) 0];
    end
    
    clear ng eg ug 
    
    wakt=wakt+toc;
    perc=j/(endd-start+1)*100;
    if wakt/60<1;  waktu=[num2str(wakt,'%.3f'),' sec'];
    else
        fwk=floor(wakt/60);
        if wakt/3600<1;    waktu=[num2str(fwk,'%.3f'),' min ',num2str(wakt-fwk*60,'%.3f'),' sec'];
        else waktu=[num2str(floor(fwk/60),'%.3f'),' hour ', num2str(fwk-floor(fwk/60)*60,'%.3f'),' min ',num2str(wakt-fwk*60,'%.3f'),' sec'];
        end    
    end    
    waitbar(perc/100,fff,['Processing sites ', num2str(perc,'%.3f'), ' % (',waktu,')']);
end
close(fff)
close all

OO=endd-start+1;   %Jumlah Titik

% ntktemp=ntk;
ntk=ntk(start:endd,1);

%Remove outside area
oaoain=1;
for oaoa=1:size(OOOO,1)
    statinsp=inpolygon(OOOO(oaoa,1),OOOO(oaoa,2),polg(:,1),polg(:,2));
    if statinsp==0
        remoa(oaoain,1)=oaoa;
        oaoain=oaoain+1;
    end
end
OOOO(remoa,:)=[];
VOO(remoa,:)=[];
ntk(remoa,:)=[];
SPCE(remoa,:)=[];
OO=OO-oaoain+1;

%Data untuk Plot GMT (Refer to ITRF)

%VERTICAL
SR0=cellstr([num2str(VOO) SPCE char(ntk)]);        
fid=fopen('V1.d','wt');
for OL=1:OO
    TOT0(OL,:)=SR0{OL};
    fprintf(fid,[TOT0(OL,:) 10 ]);
end
fclose(fid);

%HORIZONTAL
SR=cellstr([num2str(OOOO) SPCE char(ntk)]);        
fid=fopen('H1.d','wt');
for OL=1:OO
    TOT(OL,:)=SR{OL};
    fprintf(fid,[TOT(OL,:) 10 ]);
end
fclose(fid);

%Plate Velocity
%Perhitungan kecepatan Sundaland per titik
for sun=1:OO
    [DER,F]=plaref(OOOO(sun,2),OOOO(sun,1),LINP,BUJP,avel); 
    OPOP(sun,1)=DER*sind(F);
    OPOP(sun,2)=DER*cosd(F);
end

%Data untuk Plot GMT (Plate Velocity)
OOPP=[OOOO(:,1:2) OPOP zeros(OO,3)];
SR2=cellstr([num2str(OOPP) SPCE char(ntk)]);
fid=fopen('Plavel.d','wt');
for OL=1:OO
    TOT2(OL,:)=SR2{OL};
    fprintf(fid,[TOT2(OL,:) 10 ]);
end
fclose(fid);

%Data untuk Plot GMT (Refer to Plate)
OOQQ=[OOOO(:,1:2) OOOO(:,3:4)-OPOP OOOO(:,5:6) zeros(OO,1)];
SR3=cellstr([num2str(OOQQ) SPCE char(ntk)]);
fid=fopen('H2.d','wt');
for OL=1:OO
    TOT3(OL,:)=SR3{OL};
    fprintf(fid,[TOT3(OL,:) 10 ]);
end
fclose(fid);


cd(diras)

[Jeffrey Clark]

@Jenny Melynda Siagian

1. Always provide the error message(s) MATLAB provides regarding an error you want help with
2. Long programs like this should be provided as attached .m file(s)
3. Attach any input files (FIXSTA.SIG) and show any other provided inputs

[William Rose]

@Jenny Melynda Siagian,

I agree with what @Jeffrey Clark said. When I read a post with pages of code, I am inclined to move on to another post. An additional suggestion is to post the simplest possible example of what is not working as expected. You will have to get rid of all the extraneous code that does not pertain to the problem. The benefits for you are that a) you may discover the problem and fix it yourself, and b) others will be more likely to assist you.

Good luck with your work.