HP‑41C series, 78 Bytes

This program requires an HP‑41CX or an HP‑41C/CV with a time module plugged in. It works only correctly from October 15, 1582 through September 10, 4320. Ensure the display format allows for displaying 9 digits (e. g. FIX 4, the default setting). You may set flag 29 to enable the insertion of a thousands separator for every group of three digits in the integer part. Set flag 28 to use commas as thousands separators (as opposed to periods). Flag 28 and 29 set are the default setting.

01♦LBL⸆S           5 Bytes  global label requires 4 + (length of string) Bytes
02 DEG             1 Byte   ensure Degree mode is selected
─── get day number (zero-based) ────────────────────────────────────────────────
03 DATE            2 Byte   today’s date: has MM.DDYYYY or DD.MMYYYY format
04 ENTER↑          1 Byte   replicate value
05 ENTER↑          1 Byte   again because `ENTER↑` disables stack lift
   NULL            1 Byte   invisible Null byte before numbers
06 1 E6            3 Bytes  put 1,000,000 on top of stack
07 *               1 Byte   multiply so X = MMDDYYYY or DDMMYYYY
   NULL            1 Byte
08 1 E4            3 Bytes  put 1,000 on top of stack
09 MOD             1 Byte   isolate YYYY
   NULL            1 Byte
10 1 E6            3 Bytes
11 /               1 Byte   divide YYYY by 1,000,000
   NULL            1 Byte
12 1.01            4 Bytes
13 +               1 Byte   X = January 1, current year
14 X≷Y             1 Byte   swap X and Y so difference is positive
15 DDAYS           2 Bytes  difference of days – accounts for leap year
─── get fractional time of day ─────────────────────────────────────────────────
16 TIME            2 Bytes  obtain time in HH.MMSSss format
17 HR              1 Byte   convert HH.MMSSss into decimal hours format
   NULL            1 Byte
18 24              2 Bytes  24 hours
19 /               1 Byte
─── compose `day` as in specification ──────────────────────────────────────────
20 +               1 Byte   combine day number and time
   NULL            1 Byte
21 3               1 Byte   perihelion bias for zero-based day numbering
22 −               1 Byte
─── perform calculation ────────────────────────────────────────────────────────
   NULL            1 Byte
23 .9856           5 Bytes  use .985609113 for improved precision
24 *               1 Byte
25 COS             1 Byte
   NULL            1 Byte
26 −.01672         7 Bytes
27 *               1 Byte
   NULL            1 Byte
28 149,597,870.7  11 Bytes
29 +               1 Byte
30 PSE             1 Byte   delay execution for about one second
31 GTO⸆S           3 Bytes  `PSE` displays contents of X-register

NB: The actual calculator (in hardware) is too slow to perform all operations within one second. It needs roughly 3 seconds + 1 second for PSE. However, this is a circumstance outside the realm of the programming language.

C, 338

#include <stdio.h>
#include <time.h>
#include <math.h>
int main ()
{
  time_t rt;
  struct tm * ti;
  while(1) {
  time(&rt);
  ti = localtime(&rt);
  double d = 1.0 - .01672*cos(0.0174533 * .9856*((ti->tm_yday + (ti->tm_hour * 3600.0 + ti->tm_mday * 60.0 + ti->tm_sec) / 86400.0) - 4));
  printf ("%f\n", d * 149598000.0);}
}

Bash/coreutils/bc, 101 bytes

#!/bin/bash
bc -l <<<"149597870.691*(1-.01672*c((`date +%s`-`date -d 4-Jan +%s`)/5022635.5296))"
sleep .5
exec $0

This computes the offset from the 4th of January in seconds, so uses a corresponding constant to convert to radians. Half a year converts to roughly pi:

$ bc -l <<<"(365.256363/2*86400)/5022635.5296"
3.14159265361957033371

The rest of the calculation is straight from the question.

Python, 101 bytes

import time,math
a=149597870.691
while 1:print(a-a*.01672*math.cos((time.time()-345600)/5022635.53))

345600 = 4*24*3600 (four days)

5022635.53 ≌ (365.256363*24*3600)/(2π) (seconds in year/2π)

TI-BASIC, 38 bytes

Disp 25018086(59.8086-cos(5022635.4⁻¹checkTmr(83761
prgmA

For a TI-84+ series calculator. Name this prgmA. Note that this overflows the stack after a few thousand iterations; use a While 1:...:End instead if this is a problem, for two extra bytes.

This uses the perihelion on January 1, 1997, 23:16 UTC for reference, and is accurate to within a few dozen kilometers (about 7 digits of accuracy) for the next few years.

F#, 178 bytes

open System
Seq.initInfinite(fun _->
let n=DateTime.Now
(1.-0.01672*Math.Cos(0.0172*((n-DateTime.Today).TotalDays+float(n.DayOfYear-4))))*149597870.691)|>Seq.iter(printfn"%f")

This is an F# script that runs well in F# Interactive. For simplicity's sake, the "continuous output" requirement is taken to literal levels, although I did lose a byte to make the output print on a new line every iteration so that it wasn't too bad. =P

Ungolfed and explained:

Seq.initInfinite (fun _ ->            // Create an infinite sequence, with each element being defined by the following function
    let n = DateTime.Now
    let dayOffset = n.DayOfYear - 4   // Day of year returns the day as a number between 1 and 366
    let today = n - DateTime.Today    // Extract the current day, so the hours, minutes and all
    let partialDay = today.TotalDays  // Get the value of 'today' as a floating point number of days
                                      // so between 0 and 1 in this case - exactly what I needed
    // And now, the formula - note that 0.9856 has been combined with the conversion from degrees to radians, giving 0.0172
    (1. - 0.01672 * Math.Cos (0.0172 * (partialDay + float dayOffset))) * 149597870.691
)
|> Seq.iter (fun i -> printfn "%f" i) // For each of the (infinity of) numbers, print it

Python 2.4 - 158 bytes

import time,math
while 1:t=time.localtime();print(int(149597870.691*(1-.01672*math.cos(math.radians(.9856*(t[7]+(t[3]*3600+t[4]*60+t[5])/864.0/100.0-4))))))

Takes the local time and spits out the distance. time.localtime() returns a tuple and can be referenced here.

Pyth, 51 bytes

#*149597870.691-1*.01672.t*c-.dZ86400 31558149*2.nZ1

Alternate formula

d/AU = 1 - 0.01672 cos ( 2π [time since perihelion]/[orbital period] )
This formula is essentially the same as the OP's formula, except it is generalized to be able to use any perihelion as a reference date.

The OP's formula has [time since perihelion] as ( day - 4 ) and has ( 2π rad / [orbital period] ) pre-calculated as 0.9856deg/day.

In my solution I am using the perihelion closest to the Unix epoch, 2^nd January 1970.

The code

Hand-compiled to pythonic pseudocode:

#                        while 1:
  *149597870.691             print( 149597870.691 * (                 # implicit print
    -1                           1 - (
      *.01672                        0.1672 * (
        .t                               trigo(
          *                                  multiply(
            c                                    divide(
              -.dZ86400                              unixTime-86400,
              31558149                               31558149
                                                 ),
            *2.nZ                                2*pi
                                             ),
          1                                  1                        # 1 means cos
                             )))))

This is essentially just turning the following formula into code:
d = ( 1 - 0.01672 cos ( 2π (t - 86400)/31558149 ) ) * 149597870.691
where t is the Unix time.

Mathematica, 97 bytes

Dynamic[1496*^5-2501*^3Cos[.9856#&@@Now~DateDifference~{DateValue@"Year",1,4}],UpdateInterval->1]

Explanation

{DateValue@"Year",1,5} represents 5th of January this year, and ...~DateDifference~... gives the temporal distance.

Dynamic[...,UpdateInterval->1] update the expression once per second.

Java - 185 180 bytes

static void d(){while(true){System.err.println(149597870.691*(1-.01672*Math.cos(Math.toRadians(.9856*(Calendar.getInstance().get(6)+LocalTime.now().toSecondOfDay()/8.64e4-4)))));}}

This uses the fact that there are 86,400 seconds in a day and is using local time, not GMT. Output happens much more than once per second. Not sure if import statements should be included in byte count.

To include a 1 second delay adds about 26 bytes e.g.

static void d(){try{while(true){System.err.println(149597870.691*((1-.01672*Math.cos(Math.toRadians(.9856*(Calendar.getInstance().get(6)+LocalTime.now().toSecondOfDay()/8.64e4-4)))));Thread.sleep(1000L);}}catch(Exception e){}}

Java definitely isn't the most golfable language. :)

Removed a few bytes thanks to @insertusernamehere