%變量和圖形初始化，輸入氣體的摩爾數為nMoles，初始壓力為P(1)，初始容積V(1)
%氣體常數R＝8.314；給定起始總功WTotal=0；點序號iPoint=1；畫等溫線使用
%的點數NCurve=100;P-V圖第一點坐標PPlot=P(1)；VPlot=V(1)

nMoles=input('輸入nMoles: ');
P(1)=input('輸入初始壓力: ');
V(1)=input('輸入初始容積: ');
R=8.314;
iPoint=1;
WTotal=0;
NCurve=100;
%算出初始溫度
T(1)=P(1)*V(1)/(nMoles*R);                            
PPlot=P(1);
VPlot=V(1);

%在進入循環之前，設置兩個不相等的PathType和QuitType值
QuitType=4; PathType=0;
while(PathType~=QuitType)                            %在菜單上選擇“退出”之前不斷循環
    %選擇路徑類型或退出
    iPoint=iPoint+1;
    fprintf('對過程#%g\n',iPoint-1);
    PathType=menu(sprintf('過程%g:選擇下一路徑',iPoint-1),...          
        '等壓','等容','等溫','退出'); 
    switch PathType
    case 1                                          	%等壓路徑
        V(iPoint)=input('輸入新壓力：');
        P(iPoint)=P(iPoint-1);                          	%壓力不變
        T(iPoint)=P(iPoint)*V(iPoint)/(nMoles*R);        	%按新容積算出溫度
        W=P(iPoint)*(V(iPoint)-V(iPoint-1));            	%計算等壓過程所做的功
        VPlot=[VPlot,V(iPoint)];                       	%加上新的容積和壓力點，用以繪圖
        PPlot=[PPlot,P(iPoint)];
        
     case 2                                         	%等容路徑
          P(iPoint)=input('輸入新容積：');
          V(iPoint)=V(iPoint-1);                      	%容積不變
          T(iPoint)=P(iPoint)*V(iPoint)/(nMoles*R);      	%按新壓力算出溫度
          W=0;                                    	%等容路徑上所做的功為零
          VPlot=[VPlot,V(iPoint)];                    	%加上新的容積和壓力點，用以繪圖 
          PPlot=[PPlot,P(iPoint)];
          
     case 3
         V(iPoint)=input('輸入新溫度：');
         T(iPoint)=T(iPoint-1);                        	%溫度不變
         P(iPoint)=nMoles*R*T(iPoint)/V(iPoint);        	%按新容積求新壓力
         W=nMoles*R*T(iPoint)*log(V(iPoint))/V(iPoint-1);	%求所做的功         
         
        %用元素群運算求等溫路徑上的P和V，加進繪圖數據中
        VNew=linspace(V(iPoint-1),V(iPoint),NCurve);
        PNew=nMoles*R*T(iPoint)./VNew;
        VPlot=[VPlot,VNew];                           %將新的V、P點加入繪圖數據中
        PPlot=[PPlot,PNew];   
     
      otherwise
    end
    
    %畫出到目前為止的PV圖
    if(PathType~=QuitType)
        WTotal=WTotal+W;                             %將新做的功加進總功
         figure(gcf);
         plot(V,P,'o',VPlot,PPlot,'-')                     %圖形窗移前，繪圖
    end
end