I doing a three way Anova (or Ancova - I'm not sure, hence my question) in R. 
I am testing how temperature, the development stage and the size of a carcass affect the development rate of maggots.

My response variable is `Duration` (a measurement of hours) and my factors are `Size` (2 levels = small and large), and `Stage` (7 levels = eggs, 1st instar, 2nd instar, 3rd instar, postfeed, pupa and total) and `Temperature` (4 levels = 15, 20, 25, 30).

I intend to examine how `Duration` varies with `Temperature` and `Size` and `Stage`. 

I imported my data set (AnovaTWD).
One function I have tried is:

    model1 <-(aov(Duration~Stage*Temperature*Size, AnovaTWD)) 
    summary(model1)

This gave me:    
    
                            Df  Sum Sq Mean Sq F value   Pr(>F)    
    Stage                    6 7206782 1201130 149.059  < 2e-16 ***   
    Temperature              1 1924926 1924926 238.881  < 2e-16 ***   
    Size                     1   78491   78491   9.741  0.00205 **   
    Stage:Temperature        6 2500293  416716  51.714  < 2e-16 ***   
    Stage:Size               6  120539   20090   2.493  0.02367 *  
    Temperature:Size         1  140090  140090  17.385 4.43e-05 ***   
    Stage:Temperature:Size   6  184679   30780   3.820  0.00122 **    
    Residuals              214 1724431    8058                     
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Am I right in thinking that R is recognising Temperature as a continuous variable here?

Another function I have tried is:

    AnovaTWD$Temperature <- factor(AnovaTWD$Temperature)
    model2 <-(aov(Duration~Stage*Temperature*Size, AnovaTWD))
    summary(model2)  

This gave me:

                            Df  Sum Sq Mean Sq F value   Pr(>F)        
    Stage                    6 7206782 1201130 527.084  < 2e-16 ***    
    Temperature              3 2399926  799975 351.048  < 2e-16 ***    
    Size                     1  116577  116577  51.157 1.90e-11 ***    
    Stage:Temperature       18 3068185  170455  74.800  < 2e-16 ***    
    Stage:Size               6  157139   26190  11.493 6.34e-11 ***    
    Temperature:Size         3  214981   71660  31.446  < 2e-16 ***    
    Stage:Temperature:Size  18  292783   16266   7.138 1.10e-13 ***    
    Residuals              186  423861    2279                         
    ---    
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1    

I know that R is treating Temperature as a fixed factor now. However, I have been advised to treat `Temperature` as a covariate. I read "ANCOVA is easily reached using the `aov()` function using the syntax `+
variable name` to indicate that the predictor variable is a covariate."

So I tried:    

     summary(aov(Duration~Size*Stage+Temperature, AnovaTWD))
    
This gave me:

    Df  Sum Sq Mean Sq F value Pr(>F)         
    Size          1   87904   87904   4.928 0.0274 *      
    Stage         6 7223925 1203988  67.497 <2e-16 ***    
    Temperature   3 2411455  803818  45.063 <2e-16 ***    
    Size:Stage    6  143505   23917   1.341 0.2400         
    Residuals   225 4013443   17838                       
    ---    
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 


Is this now ANCOVA and treating `Temperature` as a covariate? Or is it not testing an interaction with the other factors? 

I could simply treat `Temperature` as a fixed variable as my experiment design used the 4 constant temperatures (15, 20, 25, 30). However, I know this approach is not giving all the information possible, hence I am looking into ANCOVA as it is usually used when there are both continuous and categorical predictors.